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Abstract 

This article addresses the interpolation of harmonic functions over the interior of a R'' unit sphere 
by linear combinations of fundamental-solution point-source basis functions, where all the sources are 
assumed to be outside the sphere. While it is natural to formulate approaches to harmonic approxima- 
tion, interpolation and/or boundary value problems for spherical interiors in terms of minimizing some 
standard Dirichlet integral, there is no established approach along these lines that yields interpolating 
solutions for point-source basis functions. Here it is shown that by introducing a simple weighting 
function, exact closed-form inner-products result for the expressions needed; hence, minimizing the ap- 
propriate weighted Dirichlet integral yields exact closed-form linear equation sets for the source strengths 
and the fits resulting from these equation sets match all of the prescribed values at the interpolation 
points. Further, the formalism can be extended in a natural fashion to handle interpolation fits using 
higher-order point-multipole basis functions (such as point dipoles and point quadrupoles) , so that in- 
terpolations for higher-order partials of harmonic functions can be easily implemented. Since the source 
and field points are in different domains, the fundamental-solution basis functions are bounded and can 
be regarded as defining a new type of kernel space that is related to, but distinct from, a reproducing 
kernel Hilbert space (RKIfS). This new space is labeled a Dirichlet integral dual-access collocation-kernel 
space (DIDACKS). While the corresponding DIDACKS formalisms for R^ spherical exteriors was devel- 
oped elsewhere and thus is only summarized here, a full discussion of the R^ spherical interior formalism 
is given here, including all the relevant derivations. This presentation includes numerical tests of the 
DIDACKS point source fitting method. Finally, although DIDACKS fits to data obtained from actual 
point measurements are possible, all of the required interpolation point data is assumed to be exact here 
and neither the treatment nor the consequences of measurement errors are considered in this paper, but 
ill-conditioning and condition number issues are addressed. 

Like RKHS (or least-squares collocation) fits, DIDACKS fits simultaneously satisfy two distinct norm 
minimization conditions: (1) The norm of the difference between the function to be fit and the inter- 
polating fit is minimized (the least-squares norm property). (2) Of all those functions that match the 
values of the function to be fit at the specified interpolation points, the interpolating fit is the one with 
the smallest norm (the collocation minimum norm property). The least-squares norm property follows 
directly from the DIDACKS formalism itself. A derivation of the collocation minimum norm property is 
given here, which is a generalization of the standard proofs encountered in a RKHS setting. In concert 
with these minimum norm properties, DIDACKS theory has implications for other integral and kernel 
based approaches — especially the method of fundamental solutions (MFS). Specifically, it is shown that 
natural connections between DIDACKS and MFS approaches imply that two distinct previously unrecog- 
nized minimum norm conditions hold for MFS techniques as commonly implemented; furthermore, these 
conditions imply that MFS approaches are on a more solid theoretical footing than has been commonly 
recognized. Finally, the existence of obvious connections between DIDACKS and the fast multipole 
method (FMM), as well as the boundary element method (BEM), are also pointed out, but they are not 
explored in any detail here. 
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1 Introduction 

Given a real valued function / satisfying Laplace's equation over the interior and on the boundary of an R"^ 
unit sphere, this article shows that an interpolating function Lp can be constructed from linear combinations 
of inverse-distance point-source basis functions, where the sources are located outside of the unit sphere. 

2 

Further, these interpolating solutions are obtained by directly minimizing ||/ — <y9||jj , where || • ||jj denotes 
the interior integral (II) norm — a positive definite norm defined here in a simple but nonintuitive way. It is 
shown that this norm is proportional to the Dirichlet integral of the field difference (/ — if) with an inverse 
radial distance weighting function. The developed formalism prescribes a one-to-one correspondence between 
each pair of interpolation and point-source locations so that if there are TV^ interpolation points then there 
must be iVfe point-source basis functions. (It is assumed throughout that A'^fc is a fixed finite integer.) The 
basis functions have the form \|^k^ where £fe is the distance between source point k and some field point 
that is in the interior of the sphere. These basis functions are often called "fundamental solutions" since 
they satisfy Laplace's equation for interior points, V^£^^ = 0, as well as Poisson's equation over the entire 
field and source region; moreover, other solutions can be built up from them by means of the integral form 
of Poisson's equation. Here tp :— Ik/^k, which is the discrete analog of the integral form of Poisson's 

equation, specifies the basic form of the interpolating function. The source strengths £ M are determined 
by a closed-form linear equation set that yields a ip with the desired interpolating properties. This equation 
set follows directly from the standard approximation procedure of setting the partials of 

$ := 11/ - ^llii = (/, /)„ - 2(/, ^)„ + (^, ^)jj (I) 
with respect to the source strengths equal to zero and diving by two: 



Tq = A , (2) 



which obviously minimizes <I>. Here T denotes the matrix whose elements are T^i u := {(-^'^ ^ ^fe^)ii ^'^'^ ^ 
and fc' = I, 2, 3, ... A^fc; A the vector whose elements are := (/, (■k^)^^] and finally, q the vector whose 
elements are qk ■ This interpolation procedure is made possible by the existence of closed- form inner products 
of the form ^fc and (W, (■k^)ii , where ( • , • )jj is the II inner-product associated with the II norm. 

(Subsequently, since the definition of an inner product and norm are linked in an obvious fashion, when no 
confusion can result the descriptive adjective "norm" will be used to imply both the inner-product and norm 
setting.) 

Besides the fact that the II norm yields closed form expressions for inner products of i^^ with any 
harmonic function [i.e., (fT7|) ]. it also has various other interesting properties. One is that the II norm of 
any suitable function can be expressed as the sum of a weighted Dirichlet integral term and of function 
evaluations at the origin (|I9p . while a second is that the II norm of any function can be written as the sum 
of an unweighted Dirichlet integral term and a surface integral term psp . Yet another interesting property 
emerges from the question: Does a point replication and inner product structure exist that allows one to 
reinterpret Poisson's integral formula for / in a natural way as the inner product of some fundamental 
solution kernel and / so that the end result is / evaluated at some specified point? As the derivation leading 
up to (|27p reveals, the II norm with kernel provides an affirmative answer to this question. (It is also 
worth noting that the Poisson kernel itself can be extended in a natural way that allows for a reinterpretation 
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of Poisson's integral formula as an inner-product evaluation P", p. 122]). Other interesting aspects of the II 
norm should be apparent in the sequel. 

While various procedures for harmonic interpolation have been proposed over the years that are 
largely unrelated to the approach presented here, the only direct precedent appears to be an analogous 
formalism for M.^ spherical exteriors developed by the author [19]. As one would expect, for this exterior 
problem the source points are inside the sphere and the harmonic field region is in the exterior. There 
corresponding relationships yielding closed form inner products and relationships to weighted and standard 
unweighted Dirichlet integrals are derived. It was also shown there that a similar theory can be developed 
for M"^ half-space. This exterior and half-space theory are briefly summarized in Section [21 As elaborated 
on below, the label used for this basic interpolation approach and the associated space for these and other 
geometric settings is a Dirichlet integral dual-access collocation-kernel space (DIDACKS). 

Although the goal here is to minimize the overlap between this and the companion article just mentioned, 
in order to make the present article self-contained there are three areas of modest overlap. First, given that 
there is a great deal of overall mathematical similarity between the interior and exterior M.^^ DIDACKS 
formalisms, in order to facilitate comparisons it is useful to have a ready summary of the relationships for 
both settings available and thus Section [2] fills this role. The interior relationships summarized in Section [2] 
are proved in Section[3l but the derivation of the exterior relationships are not duplicated. Second, just as in 
the companion article, it is also pointed out in Section[2]that the basic DIDACKS interpolation formalism can 
be directly generalized in an obvious way to allow for the matching of higher-order partial derivatives of /. In 
this case, the interpolating basis functions that must be used are corresponding components of higher-order 
point multipoles as summarized in Table[T] Hence from a fitting perspective, various DIDACKS multipole fits, 
including combinations of point dipoles and point quadrupoles, can easily be performed. These possibilities 
are discussed at the end of Section [51 Third and finally, just as in the companion article normalized basis 
functions are often used here, where for scalar point sources they are given by 0k /W^k^W- (Generally 

the basis functions ipk will be highly non-orthogonal.) As previously noted in the companion article, the 
use of normalized basis functions is strongly recommended for cases where large differences of physical scale 
may occur and produce unwanted effects (i.e., large condition numbers) in the solution process of ([2]). This 
occurs, for example, when scalar point sources and dipoles are both present. 

The DIDACKS interpolation procedures studied here are analogous to reproducing kernel Hilbert space 
(RKHS) interpolation techniques where an appropriate linear combination of reproducing kernel basis func- 
tions are used — Bergman kernel interpolation being one notable example. It is thus useful to briefly compare 
and contrast these two approaches. (The reader who is unfamiliar with RKHS theory may wish to consult a 
standard reference, such as [1], [16] or [17], but aside from some of the material in Appendix A these connec- 
tions are not actually an integral part of what follows and they need not be completely mastered so a detailed 
knowledge of RKHS theory is not really mandatory.) For some region C M'^, associated with each RKHS 
and its inner product form ( • , • ) is a unique symmetric reproducing kernel (SRK): K{P, Q) — K{Q, P), 
where P E M.^ and Q e M'' . Here K{X, P) satisfies the reproducing property 

(X, /) = f{P) , (3) 

where X € M"^ is a dummy argument tied to the inner-product structure. 

The inner-product evaluation property of kernels like ^ , which have one argument in one domain and 
the other argument in the complimentary domain, is called a replication property in order to distinguish it 
from the standard reproducing kernel property ([3]) , where both of the arguments involved are in the same 
domain. The term collocation property will be used to indicate either the possibility of a replication or a 
reproducing property; moreover, it will also be used in its more commonly used sense to indicate a simple 
forced matching of a linear combination of Nk basis functions to Nk specified point values. Kernels, of 
course, will be labeled by their properties so that, for example, is called a replication or collocation 
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kernel. Although it cannot be considered a reproducing kernel, since its two arguments are in different 
domains, i^^ is a bounded kernel like reproducing kernels due to the restrictions on these arguments. In 
recognition of the dual argument nature of as well as the fact that it has a replication property, it 
and like kernels will be given the generic label dual-access collocation kernel (DACK). The main focus of 
this article, of course, is the more interesting, but specialized, case where the normed space arises from a 
(weighted) Dirichlet integral (i.e., a DIDACKS), so that the kernel form employed here, , is regarded as 
a Dirichlet integral dual-access collocation kernel (DIDACK). 

Before considering the differences between DIDACK and SRK approaches, consider the similarities. Al- 
though they will not be explored in detail here, the existence of connections between DACKs and reproducing 
kernels are obvious for the DIDACK kernel form £^^. As it should be obvious from the sequel, the point is 
that the DIDACKS kernel form ^^T^ can always be transformed into a symmetric reproducing kernel form 
and then the associated replication property can be transformed into a reproducing property. For example, 
consider the II norm replication property given by ([T7)) . The DIDACK i?^^ can be transformed into a SRK 
by first using (|10p to eliminate the exterior source point in favor of a corresponding interior point and then 
substituting the result into p7|) . Dividing the resulting expression for ([17]) by Pk then yields the SRK form 
([3]) . As discussed at length in [19] , for the exterior setting a similar reformulation is possible and the result- 
ing kernel and associated weighted Dirichlet norm was first studied by Krarup [14]. The various differences 
between DIDACK and SRK based approaches were addressed in some detail in JJJ and thus need not be 
reconsidered here, but it is worth noting that vast differences are readily apparent when dipole or other 
multipole like fits are contemplated. These effects arises from the differences in action of a linear differential 
operator in the two space settings. Thus suppose that is a differential operator that is, for example, some 
combination of partials with respect to the components of Pk (or of P) , then it is easy to see that the action 
of on ([T7]) and its RKHS analog ^ are quite different and correspond to entirely different fits. Moreover, 
all of the operators considered here act with respect to source components rather than components of P^ , so 
that these differences are even more pronounced. There are also other general differences between SRKs and 
the kernels studied here, but questions of an abstract nature regarding kernel structure are not addressed in 
this article. 

Next briefly consider the general function space setting for DIDACKS. Observe that all the solutions 
to ^ are contained in the finite linear span framed by the set of basis functions used and that Cauchy 
sequences thus need not be admitted. These Cauchy sequences can be potentially problematic here; hence, 
for the main interior setting of interest the family of admissible functions is assumed to be the linear space of 
bounded harmonic functions with bounded partials in the interior of a unit sphere and on its boundary and 
Cauchy sequences of such functions are not admitted, which means an underlying pre-Hilbert space setting 
is assumed. 

Next, consider the overall structure of the paper. Earlier it was pointed out that Section [2] contains a 
summary of DIDACKS spherical exterior as well as a preview of the DIDACKS interior relationships. The 
derivations of the DIDACKS interior relationships given in Section [3] are followed up in Section [4] with a 
brief synopsis of the collocation minimum norm result for the DIDACK l/£k- This result is the analog of 
a like result from RKHS theory. A more general discussion of this result and related theorems is given in 
Appendix A. It is worth noting that the basic proof of the minimum norm collocation principle presented in 
Section[4] (which is the same line of argument that is generalized in Appendix A to include the RKHS setting) 
is very concise; moreover, the surrounding context in Section [4] and Appendix A contains various conditions 
and restrictions that should normally be stated for these RKHS proofs, but are not. Further, since there 
does not seem to be a complete discussion of all of the associated mathematical consequences and pitfalls, 
an attempt to raise most of the major points is made in Appendix A. Including the possible generalities 
and associated side issues results in a somewhat circuitous and lengthy discussion in Appendix A, but as 
mentioned in Section [4] significant theoretical and practical consequences follow from this minimum norm 
collocation principle so a thorough understanding of it is useful. In this connection, one of the primary tools 
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introduced in Section U] and Appendix A is the concept of pointwise independence and Appendix B gives a 
attempt to directly prove that harmonic functions are pointwise independent. Finally, although it is obvious 
that the minimum norm properties described in Section 3] can be used to derive various interesting energy 
norm (and weighted energy norm) inequalities, these possibilities are only briefly touched on at the end of 
Section [H 

Since the basic DIDACKS interpolation formalism developed here is new, for completeness and in order to 
give some level of assurance in this basic formalism itself numerical test results are given in Section [5l While 
all of the numerical tests in Section[5]are with respect to elementary basis functions of the form 1/ ^k, as noted 
in [19 ^ the author has had extensive favorable experience with combined point mass and point dipole fits 
(and even some with point quadrupole fits) for the exterior setting where the formalism is quite analogous; 
consequently, a duplication of these types of tests for the interior setting was not considered particularly 
relevant here. While the numerical tests were not chosen to be representative of typical applications, they 
do serve as a vehicle for the consideration of various points of interest; thus a partial list of relevant points 
suggested by the numerical results in Section O is given at the end of Section [51 As an aside, in conjunction 
with these numerical tests, a simple algorithm is described, which may be useful elsewhere, that sets up a 
semi- regular grid on a sphere (i.e., a grid where the points are more-or-less equal distances apart). 

Finally, after a little reflection, it is apparent that the DIDACKS interpolation algorithms described 
here can be used as a boundary value problem (BVP) solution technique for harmonic problems over 
not only spherical interiors, but, as elaborated on below, over general bounded connected regions that have 
suitably shaped bounding surfaces since all that is really required for the present interpolation formalism 
is that the sources be in the exterior of some sphere and that the harmonic region be contained inside the 
same sphere — in fact many of the numerical trials given in Section [5] can be reinterpreted as the method 
of fundamental solutions (MFS) results for spherical EVPs. (Using the exterior relationships summarized 
in Section [21 it is clear that a completely analogous treatment of harmonic BVP exterior problems using 
interior sources is possible.) Section [SI addresses these important MFS/BVP connections to the DIDACKS 
approach, including the significance of the above mentioned minimum norm collocation principle and the 
standard DIDACKS least-squares minimum norm principle for MFS. These minimum norm conditions imply 
that MFS approaches are on a more solid theoretical footing than has been commonly recognized and that 
this footing implies, among other things, that a greater numerical stability results than might otherwise be 
expected for deep source placements when using MFS for exterior BVP problems — which has been frequently 
observed and often commented on in the MFS literature (where it is usually noted as the confluence of good 
solutions and large condition numbers). Of course, the same source placement effect occurs for the interior 
setting and it has the same explanation. It is thus worth noting that exploring this source placement issue 
in the interior setting is one of the design criteria behind the test scenarios presented in Section [5l 

Finally, it may be useful to briefly touch on other significant topics that are not addressed elsewhere in 
the article. While there are many such topics, only three are considered below. Corresponding to each of 
these three topics is an active area of applied numerical research. As a guide to the uninitiated in any of 
these applied areas, at the end of the discussion of each of these three topics the name of journal(s) where 
relevant research articles often currently occur is given in curly braces. 

First, given the well-known historical connections of MFS approaches to the boundary element method 
(BEM) and the closely related boundary integral method (BIM) [also called boundary integral equation (BIE) 
method] techniques [5, p. 286], it would seem to be obvious that there are analogous connections between 
BEM (or BIM) and DIDACKS approaches for Laplacian problems. Moreover, for Laplacian based problems 
with BVP surfaces that have appropriate geometries, the resulting connections have significant consequences, 
including, but not limited to, implied implicit minimum norm conditions for BIM and BEM techniques that 
have previously gone unnoticed. This family of integral equation approaches is most commonly presently 
referred to as BEM approaches and they have a truly remarkable range of applicability. Some hint of these 
BEM connections is given at the beginning of Section [6] where the case of BVPs with spherical boundary 
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surfaces is discussed; however, even a preliminary discussion of BEM/DIDACKS connections is outside the 
scope of the present article. Nevertheless, before proceeding it is interesting to consider several follow-on 
questions that may, in some sense or other, prove to be suggestive: Do these DIDACKS/BEM connections 
imply that certain applications exist where a BEM-like integral approach can be developed that introduces 
a (small) layer-like region separating the sources and the boundary surface so that numerical singularities 
can be avoided? Is it possible to reformulate some BEM-like approaches to certain problems as a Dirichlet 
integral norm based approach somewhat along the lines of the DIDACKS approach (this question is raised 
by the fact that it is possible to reformulate the DIDACKS technique as a surface distribution approach). 
Above, it was noted that BEM-like approaches have a broad range of applicability — how broad is the area 
of applicability where some corresponding implicit norm criteria can be retained. Although it would seem to 
be clear that something of interest can be done for MFS or BEM problems with spherical boundaries, what 
is the most general boundary surface geometry that can be considered, which still has some relevant implicit 
norm condition? It is not clear how much can actually be done here, but these questions are raised in the 
hopes that something can indeed be done and that others will explore the possibilities that are relevant. {For 
BEM related articles see, for example, the journal Engineering Analysis with Boundary Elements.} {For MFS 
related articles see, for example, the journals Engineering Analysis with Boundary Elements, International 
Journal for Numerical Methods in Engineering or Numerical Methods in Partial Differential Equations.} 

Second, elaborating on and adapting the ideas just noted in connection with MFS/BVP approaches, it 
is readily apparent that the formalism developed here can be more or less directly adapted to problems of 
inverse source estimation for bounded regions satisfying Laplace's equation provided that the geometry of 
the bounded region under consideration is fairly benign (i.e., somewhat spherical in shape) since all that is 
really required in the present formalism is that the source region be in the exterior of some sphere and that 
the harmonic region be contained inside the same sphere. Although the region associated with the norm 
may not match the shape of this boundary, its norm minimization criteria still should generally provide 
solutions that have a significant physical grounding and it is thus not clear to what extent these shapes need 
to actually match. This issue is addressed further in Section [6l Given the general ambiguities inherent in 
inverse source estimation problems; it is, however, clear that even for spherical regions there are significant 
implementation issues and, as such, a separate article is planned to deal with them. Since inverse source 
estimation is more commonly associated with problem settings for exterior geometries, this subsequent article 
will addresses only exterior settings. Also it is clear in general that relevant points made in conjunction with 
the exterior setting also hold for the interior setting and conversely; consequently, discussions of relevant 
points will not be duplicated from one setting to the next unless clarity of exposition dictates otherwise, 
since it is assumed that these points can be carried over from one setting (or article) to another setting (or 
article). (The same comment also holds for different geometries, such as half-space settings, as well as for 
different underlying spaces.) There are nevertheless several unique problem areas of interest where sources 
are in the exterior region and harmonic interpolating solutions for spherical interiors are needed so they are 
noted here in passing for future reference. One example of where this might occur is in a physics laboratory 
or other experimental setting where a bounded static electrical (or magnetic) field test domain is of interest 
(which can be taken to be spherical or enclosed in a spherical region) and has unwanted field inhomogeneities. 
The goal, then, might be to find a set of external sources to counter these unwanted effects. Various other 
scenarios obviously exist. [Although it is outside the scope of the present article, it is worth noting that 
various ways of performing dipole I norm or II norm DIDACKS fits without assuming explicit point potential 
values exist and that one of these techniques was briefly described in [Hj.] {For inverse source estimation 
related articles see, for example, the journal Inverse Problems or the journal Geoexploration.} 

Third, a common set of problems where sources are in the exterior region and point source interpolating 
functions are desired in the interior of a bounded region are problems of a type usually addressed by the 
Fast Multipole Method (FMM). Specifically, in many engineering settings there may be numerous known 
sources in some exterior region that produce fields whose effects need to be computed in some interior region, 
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which may be taken to be the interior of a sphere (here often there are also sources in the interior region 
of interest, but their effects can easily be computed separately). The inherent difficulty is that there are 
generally so many sources in the exterior region that a direct computation of the resulting interior field is 
overly burdensome. FMMs are designed to overcome this computational difficulty and they have found broad 
usage in diverse areas such as computational field engineering and astronomy. Generally the various sets of 
sources are broken down into different regimes depending on how far away they are from the interior field 
region, so that, for example, the nearer regimes produce greater variations over the field region of interest, 
but correspond to fewer given sources. DIDACKS interpolation might profitably be used to speed up this 
computational process since a properly implemented DIDACKS fit with a reasonable number of sources 
might well be expected to replicate each of these different field variations over the entire spherical interior in 
an accurate and efficient way. To achieve this accuracy and efficiency, DIDACKS fits here clearly must be 
done for each of the needed regimes separately. Generally when FMM approaches are used the vector field 
is the quantity of interest rather than the potential itself and from Table [T] a combined scalar and dipole 
DIDACKS interpolation scheme is most appropriate for use in these cases. Thus in effect, aside from the 
additional burden of computing the potential, the only significant computational load would be the needed 
FMM computations at the few interpolation points initially required to generate the (dipole/point mass) fits 
themselves. A major point here is that there is generally a natural decomposition into near, intermediate 
and far source terms for most FMM areas of application and from the author's experience with gravity field 
modeling [12], a decomposition into these various field parts is very desirable from both a field modeling 
and computational efficiency perspective. It is clearly generally advantageous to treat parts of a field with 
different characteristics separately when using a DIDACKS based treatment. This is obviously due to very 
different spectral properties and associated attenuation attributes that each of the different parts of the field 
might have; moreover, field parts with different behaviors require quite different associated source spacing in 
order to obtain good modeling. For example, for a single source layer, positioning is important not only with 
regards to the tangential (or horizontal) source interspacing, but also with regards to the overall radial (or 
vertical) source positioning (for multiple source layers, spacing issues are even more involved, but analogous 
considerations enter). In the literature, although somewhat different proposals for simultaneously using both 
MFS and FMM approaches to related application areas have recently been proposed [HI [21], this is clearly an 
area warranting further study and in this endeavor DIDACKS theory could prove to be of general relevance. 
Although there are other pertinent points that could have been made, the above DIDACKS based approach 
just outlined was given since it is conceptually simple and may eventually prove to be of some interest 
for certain applications. {For FMM related articles see, for example. International Journal for Numerical 
Methods in Engineering.} 



2 Summary of Spherical Exterior and Interior Relationships 

This section introduces the standard notation to be used and it also outlines the basic inner-product relations 
of interest for both the M"^ spherical interior and exterior harmonic field regions. For completeness the basic 
M'^ point replicating inner-product equation needed to handle half-space interpolation and approximation is 
also presented without proof. 

Although it is slightly pedestrian, for clarity, consistency and in order to avoid various notational conflicts 
and ambiguities (such as double subscripts), in the main body of this article three dimensional Cartesian 
vectors are used and denoted by overset vectors: X — {x, y, z)"^ € (for superscript T :— transpose), while 
for n ^ 3 rt-dimensional vectors and matrices are denoted by lower and uppercase bold letters, respectively. 
(This also serves to make the 'M? nature of the relationships here transparent and serves to emphasize the 
fact that only vector analysis is required for their derivation.) For vectors, slightly different conventions are 
adopted in Appendix A as the need arises. As in |19j, the general exterior field region is: flo ■— G | 
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\X\ ^ i?o}- For the half-space case the origin is chosen in the plane dili :~ {X e R'^ | z = 0} itself and the 
field region of interest is fii :— {X G M.^ \ z ^ 0}. The general spherical interior field region is denoted by 
fi[o] := {X e I \X\ ^ i?o}, where the subscript has been surrounded by a set of square brackets that has 
been introduced here to denote the closed interior spherical region versus the subscript without brackets 
that is used to denote the corresponding exterior region. This subscript convention has been introduced 
here in analogy to the standard notation for open and closed intervals (a, b) and [a, b]: thus, for example, 
17(0) {X e I \X\ < i?o}. 

Due to the occurrence of extraneous factors of Rq, for the exterior and interior of unit spheres the 
mathematical formalism developed here assumes a slightly more elegant form that it does for spheres of 
general radius Rq. Furthermore, since it is fairly simple matter to rescale any given spherical problem to a 
corresponding unit sphere problem and then transcribe the resulting solution back, letting Rq be arbritrary 
adds little of significant value. Thus instead of dealing directly with the regions Op, 17(0) ^-i^d f7[o], the 
formalism will be presented for the following unit-sphere analogs: 



no -.^ {X e r3 I \x\^i} 
n^o] -.^{xe r3| \x\^i} 



17(0) R'l 1^1 <!} (4) 



where the commonly used hat notation for unit vectors has been mimicked. The region 17 [oj is of particular 
interest here. 

For each of these three geometries the compliment of the field region will be denoted using a superscript 
prime and locations in this complimentary region will be denoted by primed vectors. Thus a typical source 
point will be denoted = (x'j,, j/J., z^)-^ for A; = 1, 2, 3, ... , Nk, so that £k ■— \X — X'f.\. In all cases it is 
assumed that the fixed source region is bounded, closed and an interior subset of the compliment of the field 
region and this compliment is open so that l/£k is bounded. It is also assumed that X'f., ^ X'^, for all k' ^ k. 

For each of these three geometries, let f{X) and g{X) denote typical bounded harmonic functions for 
the region of interest 17, which can be either 17i, 17o or 17[o]. For the unbounded domains 17o and 17i it is 

also necessary to stipulate that / and g fall off at least as fast as 1/r as r ^ oo, where r :— \X\, so that 
\rg\ < oo and |r/| < oo as r ^ oo. Then the inner-product relations of interest for spherical geometries can 
be compactly stated in terms of weighted Dirichlet integrals over 17 by introducing an appropriate notation 
for it. Recall first that the standard Dirichlet integral over 17 is generally denoted D[/, g] — JJJ^ ^ f -^g dV. 
Generalizing this notation in a natural way to indicate not only the weighting function fi = IJ.{X), but the 
domain 17 as well gives 

D[/, <7, /i, 17] ^^^^lVf■VgdV. (5) 



Clearly D[/, 3, 1, 17] = D[/, g]. 

Consider the half-space case where 17 = 17 1 so that z ^ and < 0. Then the 17i point replication (or 
generalized collocation) property for the DIDACKS kernel I'f^^ is simply 

D[f,e-\ 1, 17i] = 2nf{x'„yl~z',) . (6) 

For the spherical exterior, where 17o denotes the field domain of interest, \X\ ^ 1 and ]X^] < 1. As an 
intermediate step it is useful to first introduce the integral (I) norm discussed in ^19j : 

(/, •= - ^// "^rirfg) da, with D, := |- , (7) 
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and where a and da have the following meaning when associated with the integral of J{X) 



II 



f{r,e, <P)da := J J [f{r, 6, 0)] ^^^sin0 ded<P (8) 

9=0 0=0 



for standard spherical coordinates r, 6 and </>, so that an evaluation on the surface of a unit sphere (i.e., at 
r = 1) is understood for the integrand of ([7]). The integral evaluation convention introduced by ([8]) is useful 
in deriving new results since, among other things, it decouples radial and angular dependencies from each 
other in surface integral expressions. 

Here the label "integral norm" derives form the fact that closed-form expressions of integrals can be 
obtained for the resulting exterior inner-product expressions that involve / and l/£fe. Specifically 

(4-1, f \ = Pkf[Pk) , (9) 

where Pk — \Pk\ with 

Pk = ■ (10) 

It also turns out that the same basic relationship between interpolation and source points holds for 
the interior case. Two other relationships exist tying the I norm to both unweighted and weighted Dirichlet 
integrals. The relationship to the standard Dirichlet integral is 

(/, g\ = (l/27r) D[/, g, 1, ho] ~ (/, g)^ , (11) 

where the second factor on the right hand side (RHS) here contains the surface inner product, which is 
defined by 

TT 27r 

(/, g)^ (l/47r) J j f{r,9, ^)g{r,0, ^)da ^ J J L=i ^ rf^d^ . (12) 

o- 6»=0 0=0 

The relationship to the Dirichlet integral with weight /io = l/*" is 

(/, g\ = (l/27r) D[/, g, /.q, ^o] = ^JJJ ^ f 9 dV . (13) 

The above relationships for Qq can be compactly stated in terms of Dirichlet and surface integrals by 
using (fTS]) to eliminate the I norm expressions: 



D[f,e-\ fio,no] = 27rPfc/(Pfe) (14) 

and 

D[/, g, 1, ho] - D[/, <7, fio, ho] + 2ir (/, g)^ (15) 

with, of course, ^0 — 1/^ and (fTO|) understood. 

For the interior of a sphere, analogous expressions exist. First, consider the actual definition of the II 
norm mentioned earlier in Section [TJ 

U.9)u--= ^fhr{rfg)da. (16) 
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While is formally similar to ([7]) they are, in fact, quite different since the domain variables for / and g 
here are understood to satisfy |X| ^ 1 rather than \X\ ^ 1. [The difference in sign between the RHSs of 
and (O obviously arises from the difference in direction of the outward pointing unit normal vectors for the 
common spherical boundary surface of the two harmonic regions.] Of course with this understanding, the 
same definition of the surface integral (|12|) and the same evaluation convention specified by ([8|) are employed. 
The same basic form of replication equation holds for the interior setting as held for the exterior, ([9]): 

{q\f\, = Pkf{Pk) , (17) 

where Pk and X'^. are again related by the same equation (fTU]) . even thought there locations have basically 
been interchanged. 

The relationship for interior regions corresponding to (llip is similar, but the surface integral enters with 
a positive sign rather than a negative one: 

(/, 5)„ = (1/2^) D[/, g, 1, f][o]] + (/, 9), ; (18) 

however, this sign difference is significant since it immediately shows that the II norm is positive definite so 
that ll/lljj = if and only if / = 0. Notice p5)) is slightly reminiscent of the first Sobolev norm where a 

volume integral of fg is added to D[/, g, 1, ^[q]]. 

The relationship of the II inner product to the weighted interior Dirichlet integral has the same weight 
function fiQ = 1/r and is also like the corresponding exterior relationship (|13p , but with an added term that 
makes it transparently positive definite: 

(/, .g)jj - (1/2^) D[/, g, ^^o, ^[o]] + fo9o^ ^JJJ r-'Vf-Vg dV +fo9o, (19) 

where fo and go denote values of / and g, respectively, at the origin: fo = /(O) and go = g{0)- For bounded 
connected regions it has long been recognized that stipulating that admissible functions vanish at some 
interior point, which is equivalent to adding a point evaluation term, transforms a Dirichlet integral into 
a positive definite norm [211 p. 109]. This clearly holds true for weighted Dirichlet integrals as well. In 
passing, notice that the singularity of fxo at the origin causes no problems since dV/r — rdrda. 

The above equations also immediately imply that compact analogs of pi)) and ([T5]) exist for the interior 
case: 

D[/, q\ MO, ^[0]] = 2^Pfe/(Pfe) - 2nPkfo (20) 

and 

D[/, 5, 1, n^o]] = D[/, 5, ^[0]] - 27r(/, g)^ + 27r fo go (21) 

where, of course, fiQ = l/r and (|10p are again understood. As an aside, observe that substituting a spherical 
harmonic expansion of / into ||/||^ shows that ||/||^ > fo and thus from (pij) that 

D[/, /, fio, n^o]] > D[/, /, 1, f^[o]] . (22) 

It is a simple matter to perform either an exterior or interior fit by substituting the appropriate inner- 
product expression [i.e., ^ or (fT7)) ] into The resulting fit can be regarded as an interpolating function 
when the Pk's are given and the X(,'s are determined by (fTO]) or as a approximation when the X^'s are given 
and the Pfc's are determined by (|10p . Additional implementation details for scalar point source fits are given 
in Section [SJ 
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Required Point Data Types 


Composite Source Types 


/ 


Scalar Point Source 


/& V/ 


Scalar Point Source & Dipole 


/, V/ & Gradient of V/ 


Scalar Point Source, Dipole & Quadrupole 



Table 1: Point Data/Source Correspondences for Spherical Geometries 



The basic formalism can also be immediately extended to interpolate for point values of the partial 
derivatives of /, which entails using other types of point source basis functions. The point is that since 
first-order partial derivatives of with respect to source coordinates and field coordinates differ by only 
a sign and differential operators with respect to source coordinates can be moved outside or inside of an 
inner-product expression at will, closed-form expressions for all of the required inner products can easily be 
obtained by taking the appropriate partials of both sides of ([9]) or pT]) with respect to source coordinates. 
The corresponding interpolating basis functions are various components of higher order point multipoles, 
such as point dipoles and point quadrupoles. There is, however, one complication: for higher-order multipole 
fits all the derivatives of the potential for all the lower orders are also required in the spherical case because 
Pk = Pk{Xl.) yields additional terms on the RHS of ^ or ([T7]) when partials with respect to the components 
of X[. are taken. This means, for example, that a dipole fit requires not only point gradient of potential 
information, but point potential information as well. For example, it is thus natural in this case to perform 
not only a point dipole fit, but a combined point mass/dipole fit. 

As a concrete example, consider how a combined point mass/dipole fit can be implemented. First, observe 
that an electrostatic or point mass dipole term is proportional to Dk-Vl'f^^ with a dipole vector source strength 
of Dk & K^. The approximating (or interpolating) potential thus becomes ip := X^^ii^fc^fc '^l ^^'^ 

the matrix T in ([2]) thus becomes a block matrix with four partitions: one corresponding to the point source 
strengths and three for the various vector components of each dipole strength. Since V^^^ = ~'^'k^k^^ where 
Vjj := {d /dx'i,, d /dy'f., d /dz'^^Y' and here components oiX'^. serve only as parameters when they occur inside 
the inner products for T and A, all X dependent derivative factors operating on i'^^ that occur inside inner 
products can be replaced with derivative factors; consequently, these differential operators can be moved 
inside or outside the inner products entirely as desired so that all needed inner products for T^^k' and Ak 
can be easily evaluated in closed form. The associated point measurement quantities and point sources up 
through composite quadrupole fits are displayed for both the spherical interior and exterior cases in Table [T] 
Analogous possibilities exist for half-space fits, where matters are even simpler since in this case for a dipole 
fit, quadrupole fit or other higher-order fit (or any combination thereof) there is no need to include the lower 
order point source types. 



3 Spherical Interior (il[o]) Relationships 

This section derives the three primary relationships for the II norm: (fT7|) . ([T8|) and (fT9|) . Although they are 
of secondary interest here, the corresponding relationships for the I norm {([9]), pT|) and (|13p} can be derived 
in a similar fashion. In what follows liberal use will be made of the evaluation convention implied by ([5]) and 
all functions considered will be assumed to be harmonic in the domain f2[Q] and have bounded partials. 
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First consider the derivation of ^T7\ . From definition (fTOl) it follows immediately that 

(5- /)„ - i^JJf9da+^JJ rg'D.U) da+^jj rfDAg) da . (23) 

Since V^/ — and V^(? — 0, applying Green's second idcntit}0 with dS — da and d /dn :— "Dr shows 

gVrif) da = ff fVrig) da (24) 



and thus that can be rewritten as 

(.9, /)„ = ^// /{2 2)'-('^5) - .9} da . (25) 

(T 

The strategy will be to show that when g is replaced by l'^^ in (pS)) . then ([T7)) follows from Poisson's 

integral formula for the interior of a sphere [13j p. 526]. Towards that end let := | and cosipk ■= 

{X ■ X'j,) / {r r',.) . Then substituting = l/^^r^ + rj,^ — 2rrJ, cos 1/;^ for ,g into (^5]) and performing the 
indicated evaluation leads to 



f)^^= l-JJ ^-!^L_J11 da, where Lk:= sJl + r'^'~2r[cos4'k. (26) 

(7 

From (fTU|l — so that can be easily reexpressed in terms of Pj^: 

K\ /)„ ^ ft (1^)// ^ P') 

From Poisson's integral formula for the interior of a sphere, the RHS of (P7)) can immediately be recognized 
as Pk times f{Pk), which proves ([17]). Although < cx) [and thus from (fTO]) that Pfe > 0] has always been 
assumed, for completeness the relationship corresponding to p?|) for the case Pk = will be derived at the 
end of the section. 

Next consider the derivation of (fT8|) . From Green's first identitjH 



Vf-Vg dV = 1 1 gVr-f da . (28) 

"[0] 

Thus (dH) follows directly from (^5)) . and 

Proving is a bit more involved. First recall the standard identity 

VV'/> = ^V^^ + (/)VV + 2VV' -V^ . (29) 

Setting ^ = fg and cj) = r^^ in ([29]) gives 

V'(/5/r) = /ffVV-i + r-^VVs + 2 Vr-^ • V/5 . (30) 



' llLi'^^'^ + ■ V^) dV = Jfg^ <i>^ dS 
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Using dH]) to reexpress V^fg as 2 V/ -Vy , along with the facts that Vr^^ • V = -{X/r^) ■ V = -{l/r^)d/dr 
and that V^(r^^) = —4:Tt6{X), where 6 is the Dirac delta function, allows ([50]) to be rewritten as 

-V/.V,= V(/^)+2./,^(X) + l^ (31) 
r 2 \ r J r'^ or 

Next consider Gauss's divergence theorenjf] for the region of interest with F = V(/5/r): 

i2( f9\ ff^ffa 



// a^l^f M'^ (32) 



"[0] 

Integrating ((3T|) over 17 [q] and using ([32]) gives 



1 



^„,,,_i rra „ , rr [ d{fg) 



r-^Vf-VgdV^-JJ -^fjda + 2n^g^ + JJ ^ J ^dr^da, (33) 

where dV^r^drdcr was used on the last term on the RHS here. [With respect to the last term on the 
RHS of ((33)) . since the r limits for the integrand occuring inside the integral with respect to dcr are fixed it 
is clear that the evaluation convention implied by ([5]) is immaterial for this particular integrand.] 
Since T>r{fg/r) = —r^^fg + r^^Dr{,fg) , ([55]) can be rewritten as 

D[/, g, 1, 17[o]] ^-lJJf9d<T+lll da + 2nfogo + JJ fg da ~ Anfogo 

I fr " 

= 2 JJ '^rirfg) da - 2nfogo , 

a 

which immediately shows ([T9|) from the definition of the II norm, (fT6| . 

Finally consider the analog of (fT7|) for the case Pu = 0. First obverse that since X'y. and are related 
by (Uni) 

lim = lim — = 1 , (35) 

P.-o Pdk P.^o \Pf^X^PkXk\ 

thus allowing interpolation at P^. = corresponds to adding a constant basis function to tp. Thus using p6p . 
consider the following inner product 

a^f)u= i^JJ^rirf)da= l.JJfda+l-JJl)Af)da. (36) 

Next consider a spherical harmonic expansion of /. From the spherical harmonic orthogonality conditions 
it is obvious that 

^JJfda=U, 

a 

while from Green's first identity or Gauss's divergence theorem it is clear that 

l-[[vAf)da^O. 
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Hence the inner product analog of (fT7|) for the case Pk = is 

(1, /)„ = /o, (37) 
which also holds for the case / = 1, as a direct substitution into shows: 

(1, 1)„ = 1 ■ 

(If the evaluation of any inner product expression is ever in question, then the definition of the inner product 
itself should be used for the evaluation.) 



4 Synopsis of the Collocation Minimum Norm Property 

This section gives an overview of the collocation minimum norm property that includes an indication of its 
utility. A fairly detailed mathematical discussion of the collocation minimum norm property for a general 
setting with various proofs and counter-examples is given in Appendix A; however, this appendix is intended 
for a limited audience. The discussion here and in Appendix A has implications for geophysical collocation 
(GC) (perhaps more commonly called least-squares collocation [U]), general kernel fitting theory (specifically, 
RKHS collocation fits), interpretation of inverse source theory and MFS like approaches (including BEM 
and related approaches for certain relevant geometries) . Special care and attention is needed here since this 
result must be extended to DIDACKS theory, but the generalizations of it considered here may also prove 
to be useful in other contexts. 

For the II norm, the collocation minimum norm property for scalar point sources over fi^o] can be stated 
as the fact that of all admissible functions g such that 

giPk) = Ck for fc = 1, 2, 3, . . . , iVfc , (38) 

where the Ck are arbritrary bounded constants, a DIDACKS scalar point source fit with coefficients qk 
specified by ([2]) is the function of smallest norm: \\g\\-^j > \\<p\\-^y Obviously this means that there are two 
distinct minimum norm conditions simultaneously satisfied by a DIDACKS fit that is based on l/ik' 

• 11/ — (y^ll is minimized directly by ([2]). 

• \\ip\\ is minimized by (j38p . 

(Trivially, if > then \a\ > \b\ and conversely, so that "yslP and \\(p\\'^ are also concurrently minimized.) 
The first of these two norm properties is obviously naturally satisfied by any DIDACKS fit. As shown in 
Appendix A, the second minimum norm property also holds simultaneously for DIDACKS scalar point source 
fits over exterior geometries. RKHS fits that meet certain minimum requirements also display both of these 
two minimum norm properties. For a general RKHS setting, although they arc usually not addressed in 
the general literature in conjunction with demonstrations of the collocation minimum norm property, the 
precise theorem requirements necessary are explicitly spelled out and discussed at length in Appendix A and 
the topic of counter-examples is also raised there in order to bolster the case for normally including these 
conditions. (Unfortunately no references known to the author discuss these conditions.) 

There are two primary lines of argument associated with the collocation minimum norm property (|38p 
that are of interest here: 

• Type I proofs, which are based primarily on various properties of the fit ip itself and 

• Type II proofs, which are based primarily on the properties of the function to be fit, /, 

as well as its associated function space. 
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Appendix A details various aspects of each of these approaches in turn. In order to clarify the various 
relevant aspects, it is useful to give a brief overview of each of these overall lines of development. Type I 
arguments are generally used to prove RKHS collocation minimum norm property (i.e., GC) analogs of (|38p 
[17j . Type II arguments are generally found in conjunction with proofs of associated results that are often 
found in mathematics monographs and texts and, in a sense, are otherwise hard to characterize. Both of 
these two lines of argument have been slightly abused in much of the literature, but the first one perhaps 
much more so than the second one. The main abuse of for both lines of argument in the existing literature 
has been to state results as though they were perfectly general without warning the reader that, while they 
may be true in the particular context under consideration, they cannot be applied to other contexts without 
closer examination and may, in fact, well not hold at all in these other contexts. For example in the main- 
stream GC literature, historically in proofs of collocation minimum norm results it has been the convention 
to implicitly assume that the analog of exists when, in fact, as it should be apparent from the sequel, 
this actually deserves some sort of proof and, furthermore, it clearly is not true for some common RKHS 
settings of interest — see, for example, Moritz's Advanced Physical Geodesy, which is a well written and well 
respected GC reference [ITl pp. 207-220]. (Helmut Moritz, along with Torbin Krarup, actually founded GC, 
hence Moritz is rightly considered the leading expert in the field.) Generally the places where the existence 
of enters as an unstated assumption in these sorts of proofs is not obvious since the proofs themselves 
are usually longer or more abstract than they need to be and, as the line of proof given immediately below 
shows, the implicit point of entry of this assumption itself is not at all obvious — except after the fact and 
when one knows to look for it. 

Before proceeding, consider the overall way that the above two types of arguments will be employed. 
The overall distinction between the two types of arguments above may be slightly unclear for now, but this 
distinction should be clearer from the sequel since these arguments are used in a fairly clearcut way. First, 
a Type I argument will be used to show that when exists the collocation minimum norm property 
exemplified by ([55)) holds for DIDACKS scalar point source fits and that, moreover, from results in the 
literature can be shown to always exist for this case. All this results in Theorem[Tl which is generalized 
in Appendix A. In general, the major problem here is that due to properties of (or restrictions on) the 
space of admissible functions ([38|l may represent an overdetermined or nearly overdetermined system and 
the easiest way to prove exists is by examining the properties of the space of admissible functions itself. 
This Type II line of thought culminates in Theorem O Since, as discussed below, scalar point sources are 
linearly independent and since DIDACKS point replication conditions exist, in Appendix A it is shown that 
M.^ harmonic functions over Hq, which is the main geometric backdrop for GC, satisfy a certain property 
(namely pointwise independence), which, in turn, implies that if no factors including differential operators 
are involved, then the GC analog of always exists. Thus, at least for this special case, one can then say 
that the usual GC collocation minimum norm proofs do indeed hold. 

Thus first consider an outline of how (|38p might be demonstrated using the Type I approach. For 
convenience it will be assumed that Pk =/= 0, although it should be apparent that the same line of argument 
holds when one of the points in (l38|) is zero. Here, of course, it is assumed that / is a fixed specified function 
and that ip = J2k=i Ikl^k is also a fixed specified function, where the source strengths qk are determined by 
([2|). Hence ([2]) can be written out explicitly as 

T.^^'k\tl')n1^'-{^l\f)n^ (39) 

k' = l 

or as 

{(■k\ V)u = i^k\ /)„ , (40) 
since the qus are assumed to be fixed by ^ itself. Because Pk> (HZ]) then immediately implies 

^(Pfe) = /(Pfe) . (41) 
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Let Ck '■— f{Pk), then (|¥T|) directly shows that ip satisfies the point matching requirement of ([55)1 . 

It remains to be shown that no other function that satisfies the same matching condition has a larger 
norm than ip. First, for convenience rewrite ([T]) as 

m := 11/ - ^llfi = - 2(/, + Ml > . (42) 
Multiplying (|40)) by qk and summing over k yields 

Ml = /)ii • (43) 

2 

Using 2{(p, = 2||iy9||jj in produces 

*(/):= 11/ -^Ilii = ll/Ilii-Il¥'llii>0 or |1/|1„ > . (44) 

Next let g be some arbritrary admissible function that satisfies the point matching equation ([38]) at 
the specified locations. Consider a DIDACKS fit to g with basis functions and coefficients g^: p' := 
Sfc^i I'k/^k- Since the fit of to g is completely analogous to that of pi to /, it follows that a relationship 
corresponding to ([33|) also obviously holds for g and p?': 

Hg) ■■= 11.9 - ^'llii = ILgllii - ll^'ll'i > O or ||.g||„ > ||^'||„ . (45) 
Next consider the actual form that the fit p' takes. For fc = 1, 2, 3, • • • Nk, since 

g{Pk) = Ck = f{Pk) 

it follows from P7)l that 

9)u = (^k\ /)„ 

and thus it follows from p9|) that = qk so that p>' = p). (|45| then shows \\g\\^^ > \\p\\-^^, as desired. This 
proof is a streamlined and much more transparent version of a standard type of proof of the collocation 
minimum norm property that is found, for example, in discussions of GC. 

Consider a general RKHS setting. On first examination, it might seem that this line of proof just 
given is perfectly general and can be directly adapted to any RKHS setting without the need for additional 
stipulations or conditions. Upon closer examination and upon considering possible counter-examples — as 
in Appendix A — it is apparent that at several places T"-'^ was assumed to exist and that this is needed 
as a separate assumption. As discussed in Appendix A an independent proof exists in the literature that 
scalar point source basis functions are independent, which implies here that T^^ always exists, so this is not 
needed as a separate assumption for the special case dealt with above. (Here linear independence of a set 
of functions {/j(^)}^i means that X^jLi ^jSj{^) = implies Cj = for all j and all X The relevant 
points pertaining to linear independence of functions can be gathered, for example, from Korn & Korn at 
various places p. 488, p. 253].) Since it is obvious that the above line of argument also holds for any 
appropriate DIDACKS setting, this can all be summed up in the following theorem: 

Theorem 1. A finite R'^ DIDACKS scalar point source (i.e., point charge or point mass) fit p exactly 
replicates the data values that result from the point evaluation of some specified function to be fit f ; moreover, 
the resulting fit minimizes not only \\f — p\\, but also \\p\\ is less than or equal to the norm of that of any 
other admissible function that also matches the specified data point values, where \\ ■ \\ is an appropriate 
DIDACKS norm. 
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In general it is obvious that if the dimension of the family of admissible functions is N]^ or less, then 
the point matching condition (|38p determines (at most) one unique function so that at best (f = g — f and 
thus $(/) = 0. When there is only one unique function specified by (|38p the two minimum norm conditions 
become the trivial statements that ||/ — (^|| = ||0|| = and \\g\\ ~ W'fiW- It is useful for having a way of 
preventing this, but the concept of linear independence of functions is not strong enough. Thus a new 
condition is needed here to quantify the relevant linear-independence like property of the class of admissible 
functions. Thus, consider the following definition that is stated for convenience in terms of R'^ points, but 
holds for scalar valued functions over R": 

Definition |4ll A class of scalar valued functions 3^ is said to be uniformly N pointwise independent over 
a region if for iV + 1 arbritrary points Xj e fl and arbritrary bounded constants Cj , there always exists 
a function f 3^ such that f{Xj) = Cj for j — 1, 2, 3, . . . , iV + 1. For any finite N, if 3^ is uniformly N 
pointwise independent, then 3^ will be called simply uniformly pointwise independent. 

The relationship of this definition to linear independence of a set of functions will be briefly considered in the 
next paragraph after a few other points are made here. As a slight abuse of nomenclature, if no confusion 
can arise a function / may be labeled uniformly (TV) pointwise independent rather than the class of functions 
that / belongs to. Clearly, if the class of admissible functions is uniformly N pointwise independent, then 
it has a dimension of at least N and there is more than one function in the class of admissible functions 
that satisfies (|38|) when Nk = N. Finally as an aside, although it might seem obvious that most classes of 
functions satisfy Definition 3)1 in one sense or another so that it is a trivial property, it is actually a fairly 
strong property. Suppose, for example, that a family of functions is uniformly pointwise independent over f2 
and that a complete orthonormal set of basis functions exists for this family over the same region, then it is 
obvious that at any arbritrary point of all of these basis functions cannot simultaneously be zero (i.e., at 
any point at least one of the basis functions must be non-zero). It is also easy to show that a restricted class 
of functions that necessarily is either even or odd with respect to any of its variables cannot be uniformly 
pointwise independent. Also the demonstration that a class of functions is uniformly pointwise independent 
usually requires a constructive proof. 

The concept of a uniformly N pointwise independent span of functions is stronger than that of set of span 
of A^-|- 1 linearly independent functions. One type of counter-example will have to suffice here. The functions 
under consideration need not be harmonic, thus let fj € {fj{X)}^Ji^ be a function that has the value 1 at 

Qj and that rapidly falls off to zero away from Qj so that fj{Qj') — for j' ^ j. This set of functions is 
obviously linearly independent. However, for most of these types of functions one would generally expect 
to find some points in f2 such that fj{X) = for all j and when this happens the set 3 := linear span of 
{fj{X)}^J^-^ is obviously not uniformly A'^ pointwise independent. 

Next briefly consider Type (II) approaches to collocation minimum norm related results. One common 
feature of these approaches is that while mathematicians generally use them correctly, the results obtained 
are, in the end, only loosely related to the least-squares collocation minimum norm principle. The Type (I) 
approach used above to obtain the collocation minimum norm result given in Theorem [T] will be the main 
one used here and in Appendix A, so the Type (II) approach used below will be primarily employed as a 
tool for obtaining a different perspective. This idea is that if the class of admissible functions is uniformly 
pointwise independent (or uniformly N pointwise independent) then the existence of a point replication or 
reproducing kernel has strong consequences, that include, among other things, the linear independence of 
associated basis functions. Given, then, that the basis functions are necessarily independent the above proof 
suffices to show \\g\\ > \\ip\\; hence, the issue of linear independence of basis functions is the only one that 
will be followed up on. 

Since little additional effort is required it is useful to consider the case of both point replicating and 
reproducing kernels simultaneously over some domain il; hence, consider a linear expansion in terms of 
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basis functions of two arguments: ip :— 1kB{X, Pk), where X e is regarded a an independent 

vector variable and Pfe g is regarded merely as a vector parameter. Here (point) replication conditions are 
generalized so as to handle both point replicating and reproducing kernels by assuming that for B = B{X, P) 

{B, /) = h{P)f{P) (46) 

where h{P) ^ — or more specifically for cases of interest that h{P) has the form h{P) = A or h{P) = P for 
some nonzero A. Thus, for example, (|17p can be reexpressed in this form by first using (jlOp to eliminate the 
source location X' in £~^, which results in the required B form with h{P) := P. [Here the case P = can be 
handled by using ((37|) .] Likewise, for a reproducing kernel K{X, P), since {K, /) = f{P) the identifications 
B{X, P) := K{X, P) and h{P) = 1 suffice. 

In Appendix A the basis functions are defined in a slightly more general way where the vector parameter (s) 
need not be explicitly introduced as an argument. Thus, in this sense, let Bk{X) := B{X, Pk) where the 
points Pfe are assumed to be distinct, then the following relevant theorem holds: 

Theorem 2. // the family of admissible functions is uniformly pointwise independent and if there is a point 
replicating or reproducing kernel characterized by [J^ l then the set of basis functions {^fcj^i are linearly 
independent. 

Proof. Suppose to the contrary that the Bk's are linearly dependent and thus that some set of q^s exist 
such that X^fe^i Bk{X) qk ~ 0, where at least one of the qk^s are non-zero. Specifically, let qk' ^ and then 
reorder the basis functions so that Bk' is the first one and thus qi =/= 0. ((46|) implies that '^k=i IkiBk, f) = 
Qkh{Pk)f{Pk) = 0. Since / is uniformly pointwise independent, choose it such that f{Pk) = 1 if A; = 1 
and f{Pk) = 0, if k > 1, then X^aS^i 1kh{Pk)f{Pk) — <lih{Pi) ^ 0, which is a contradiction. □ 

Finally briefiy consider the significance of the collocation minimum norm property. As noted above, 
the collocation minimum norm property used in conjunction with various covariance norms over f2o is a 
cornerstone of GC [17]. It is also well known result from general RKHS theory within various other contexts, 
such as analytic function theory and functional analysis: given certain conditions, of all those functions that 
match some specified function (/) at Nk points, a linear combination of Nk (symmetric) reproducing kernels 
is the one that has the smallest norm. There are also two broad areas where this minimum norm property 
strongly suggests that DIDACKS theory has results of direct interest for other mathematical specialties: 

• As previously mentioned, for MFS and BEM-like methods DIDACKS theory implies that for many 
Laplacian problems there is a not only an implied minimization of ^ , but also that the MFS or BEM 
solution, ip, is the one of smallest possible norm; i.e., any other harmonic function g that matches the 
prescribed boundary conditions (which are generally discretized or otherwise incomplete) must satisfy 
\\g\\ > ll'/'ll- Observe that this condition provides a very strong interpretational basis to MFS and BEM 
methods for appropriate geometries that, so far, has gone unnoticed. Moreover, it probably implies 
that there is a certain inherent numerical robustness in these approaches. Furthermore, one could 
argue that even when these approaches are based on some other criteria than the direct matching of / 
at prescribed points, the possibility of such a match is always implicitly there and thus the associated 
theoretical implications are also, in some sense or other, present. 

• It seems likely that various (weighted) energy norm inequalities can be derived and that they have 
both theoretical and practical importance (c.f., Raleigh-Ritz approaches). 
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Connections of these minimum norm conditions to MFS approaches will be dealt with briefly in Section [S] 
and those of BEM-like approaches are presumably analogous, which leaves only the energy inequalities to 
discuss. 

It is clear that DIDACKS theory is useful for deriving various energy (norm) based inequalities from an 
examination of the interactions of ((T4| . (fTSl) . pO]) . pT|) and p2|) in concert with these two minimum 
norm results. A single example here will have to sufflce. Thus consider the exterior of a sphere (i.e., the 1 
norm over D,q) and a harmonic function / that falls off sufficiently fast as r — )■ oo. First, from the collocation 
minimum norm property applied to the exterior 1-norm setting 



where q and (p^ are assumed to be fixed by 12]) for the 1 norm and thus, for example, from ([9]) A is the vector 

whose components are Ak = (i^^ , f)i = Pkf{Pk)- Since here D[/, /] = D[/, /, 1, > D[/, /, /iq, ^o] , from 
inequality ((47|) it follows immediately that 



For Nk less than or equal to three can be evaluated in closed form and thus a closed form expression 
for inequality (j48p can easily be found, while for larger Nk a numerical lower bound can be easily computed. 
(The explicit determination of these 1-point, 2-point and 3-point forms is left to the reader.) 

5 Implementation Issues and Numerical Tests 

This section presents and draws conclusions from numerical interpolation examples based on ([2|). While 
the interpolation region of interest is always taken here to be the entire interior of a unit sphere, it is 
also obviously possible with the present approach to interpolate over only some portion of this interior. 
As mentioned in Section [U only interpolation using scalar point sources will be considered here since the 
author has had extensive experience with dipoles and quadrupoles for the exterior of spheres and for 

half-space ^lOj and this experience strongly indicates that with a little care others should be able to 
apply DIDACKS theory to standard LLSQ dipole interior problems without untoward difficulty; moreover, 
the difficulties that might normally be expected arc mostly the same ones encountered in performing LLSQ 
scalar point source DIDACKS fits. 

Before detailing the actual specifications of the test functions to be fit (/) and of the interpolating 
functions ((/?) used to match them, several computational and conceptual details need to be briefly addressed 
since the matrix T in ([21) often has large condition numbers. After this discussion various other side issues 
will also be addressed before the actual numerical examples are considered. 



Large condition numbers are not unexpected since point source problems are known for their ill-conditioned 
nature; moreover, the confluence of large condition numbers and good results has been long noted in conjunc- 
tion with the somewhat analogous MFS approach [3]. The magnitude of this confluence and the necessity 
of pairing a suitable f to /, based on the properties of /, are two of the major results that follow from the 
numerical examples. A simple explanation for this confluence exists and is given next. Here the condition 



||/i|?>||^Ji; = q^Tq = q^A 



(47) 




(48) 



Side Issues 
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number {C^) is taken to be the largest eigenvalue divided by the smallest eigenvalue and absolute values are 
not needed since all of the eigenvalues encountered here are positive. 

First, the same basic confluence of large condition numbers and good fits occurs for point source fits over 
any region and since physical geodesy and MFS approaches commonly use exterior spherical settings and 
the language is a little less cumbersome for these settings, for simplicity VLq will be assumed in the present 
paragraph. (Recall that in this setting X > 1 > X^. > 0.) Second, it is easy to see that it is possible for 

if to match / even in the presence of large C^/^'s since ([2|) directly minimizes ||/ — ip\\-^ and there are no 
constraints put on the source parameters themselves. Hence, the values of \qk\ may be large and yet the qk 
may still produce a (f that matches / closely. Third, when normalized basis functions are used, the resulting 
T matrix will be denoted T and the mathematical form of its elements are suggestive: 

n,,:=ii-\i-}\/{\\£-X\\i-,X}. (49) 

From an examination of the explicit closed-form expression for (j49p [which can be easily determined from 
(|13p and (jl4p ] it is obvious that deep source placements correspond to large inner products when k' ^ fc, 
which translates directly into large C^'s. Fourth and finally, if / is smooth, then the basis functions 
that are used to match / must also be smooth, which means that the sources must be deep. One important 
point here is that if / is not smooth then shallow point source placements may be called for and, in this 
case, the C#'s may be much smaller. 

Consider this last point. In the MFS literature it has often been observed that the deeper the sources are 
the better the results and it has occasionally been speculated that this is a general characteristic. To see that 
this is not so and that it is dependent on the smoothness of /, only one simple counter-example needs to be 
considered. Thus consider a function / to be fit with the form / = qf/\X — X'j:\, where g/ is some appropriate 

constant and X'j: is fixed in what follows. Suppose, that (p has the usual form: ip := X^fc^i Qk/^k and that 
the "depth" of the sources is adjusted simply by simultaneously rescaling the ^^'s by Xj. ^^k^ so that 
the various source locations determine fixed directions even as A is varied. If for some fixed constant and 

some fc', X'jjr'j = X'^,,/r'^,, then as A is adjusted the miss-match in the fit will improve or worsen. In fact, 
11/ — (^11 J ^ as rj. — > r'j either from above or bellow and, to within allowed numerical accuracy, a perfect 
fit results when rj. — t'j- A similar example obviously results when / is defined by an appropriate surface 
distribution of points. Observe that when expanded in terms of a spherical harmonic series qf /\X ~ X'j\ has 

spherical harmonics of all degrees and orders and that as \X'j\ 1 the relative size of the higher degree and 
order terms increases. In short, to properly fit a choppy / requires many sources on the surface, but to fit 
a smooth / requires fewer, but deeper sources. In order to underscore and explore the phenomena of deep 
source placement here smooth test functions for / are used. While the above discussion may seem obvious 
after a little reflection, what is unexpected is the actual magnitude of the C#'s that are encountered here in 
the presence of good fits. 

Several questions naturally arise here: What is a precise quantitative statement of the underlying phe- 
nomena, which has been noted in the MFS literature, that clearly links deep fundamental solution source 
placements, large C#'s and smooth functions to be fit? Moreover what are the actual limits, in both a 
practical and theoretical sense, of actual depth of point source placements here? For the particular examples 
considered, the issue of deep source placement is explored to the numerical limits allowed by the hardware, 
software and operating system used in this study; however, it would seem quite likely that other operating 
environments allowing for more significant digits and thus deeper source placements might yield even better 
results for the test cases studied here. [Here, it seems reasonable to speculate from the source statistics 
given in Table IH that some sort of aggregate multipoles are being formed in the usual sense of the standard 
limiting process employed by physicist to define specific multipole types (for an electrostatic dipole with a 
specific orientation, this process consists of letting two point charges of equal and opposite sign approach 
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each other along a specific direction, while simultaneously letting the magnitude of the charge approach 
infinity in such a way as to give a well defined dipole limit). Following up on these possibilities is beyond the 
scope of the present article and it probably will not be fully addressed by the author in the future either.] 

One more surprising aspect associated with large C^'s uncovered by this study is worth explicitly noting: 
After the onset of numerical roundoff due to large C^'s, the actual modeling error of the resulting fits 
degrades gracefully. [This is undoubtedly due not only to fitting properties of the II norm, but also to 
properties of the Householder triangulation algorithm used in solving ^ and it strongly hints at underlying 
numerical robustness.] 

Obviously, to obtain reasonable results in the presence of large C^^'s attention must be paid to equation 
solving techniques employed and to computational word length(s) used. Toward this end all computations 
were performed using a 128-bit word length (33-34 significant digits). (On many systems 128-bit words have 
fewer significant digits than this.) The linear equations themselves were then solved using Householder trian- 
gulation. The eigenvalues, and thus the C^'s themselves were computed using singular value decomposition 
(SVD) software implemented with the same word length. 

Next consider the overall strategy employed for placing the interpolation points. Observe from the 
collocation minimum norm result for scalar point sources — Theorem [1] — it is obvious for ^Iq that when the 
interpolation points Pk are placed on a sphere of radius Rp < 1 that, for well-behaved functions /, ip will 
match / on and over the interior of this sphere of radius Rp if the interpolation points are sufficiently dense 
on this sphere. As an aside, the properties of / that are necessary to insure that it can be considered 
well behaved can be made fairly precise. From subsequent considerations based on the Nyquist frequency, 
it is apparent that if / is to be matched over some grid with a specified angular spacing, then it must 
have no frequency content above a certain limit (i.e., its spherical harmonic expansion must not have terms 
larger than some specified degree and order). Conversely, if / doesn't have such frequency limits, then 
from Theorem [S] in Appendix A and DIDACKS theory examples displaying arbitrarily large modeling errors 
can be constructed. Also by harmonic outward continuation it seems reasonable to expect that (p and well 
behaved f's will match to a certain extent over the region characterized by 1 > r > Rp. In fact, contingent 
on a good match for the region 1 > r > Rp, this interior spherical surface interpolation point placement 
strategy is perfectly reasonable in practice. (At least for the functions tested here, this is indeed found to be 
the case.) While this strategy will clearly be generally suboptimal in the sense that for any specific / under 
consideration if, as discussed immediately below, a nonlinear least squares (NLLSQ) fit of f to f is performed 
a better fit can almost always be found; there are many circumstances where a generalized LLSQ fit based 
on ([2]) is not only adequate, but preferable. Not only are these NLLSQ fits much more difficult to handle, 
but they introduce a number of unique pitfalls of their own due to the real possibility of large condition 
numbers and their ramifications, which include the implications of extreme nonlinearity. One major point 
is that, when the spectral properties of / are well known, an appropriate fixed grid can be selected and the 
resulting fitting properties of ^ can be exhaustively analyzed so that a repeatable fitting algorithm can 
be produced, which is only possible for NLLSQ fits if the very greatest of care is taken. (Observe that the 
counter examples discussed above based on f = qf /\X — Xj:\ do not apply in cases where spectral limits are 

imposed since, for finite nonzero the form 1/\X — has an infinite spherical harmonic expansion.) A 
generalized LLSQ fit based on ^ also has the following advantages: 

• A natural correspondence to spherical BVPs is set up. 

• Interpolation depth (Rp) and angular spacing parameters can be easily understood and implemented. 

• Regular samplings or griddings of point masses have a long history with a rich associated (geophysical) 
literature, which aside from elementary considerations is not required here, but that can be called 
on as needed. (A reference to this literature can be found in |19|.) An underlying assumption of 
this literature seems to be a completely regular gridding and, as such, the theory is unambitious for a 
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plane, but the full implications for a spherical gridding are unclear. The development here is completely 
independent of the ideas in this literature and, as such, it would be interesting to consider the of ideas 
of this literature within a DIDACKS context, but this can be said of many other research areas as well. 

• The types of functions / that can be fit for a given angular spacing can be easily characterized and 
understood, so that an appropriate point source spacing can be chosen for any specified /, or if / 
cannot be fit, then the reasons for this can be easily understood. 

A NLLSQ fit of to / entails iteratively adjusting the source locations to minimize so that the resulting 
Pk 's will end up on a spherical surface for only very specialized types of /'s. A further consideration of NLLSQ 
fits, however, makes it obvious that there are no universally valid interpolation point placements, which are 
optimal for all /'s, since for any initial Pk values an / can be found that will produce a lower value of <& 
when the Pk's are changed. To see that this is so, suppose that the Nk fixed locations Pk have been chosen 
as a hypothesized good universal source placement candidate, so that they also serve as an initial NLLSQ 
guess. Further suppose that / itself is specified by Nk sources, which are all at distinct locations form the 
chosen values for the Pk- In this case an adequate NLLSQ fit will obviously produce a considerably better 
value of If than a simple linear fit will. This might seem to argue for the necessity of performing NLLSQ fits 
and, in fact the spherical exterior formalism was developed by the author over 25 years ago expressly for the 
purpose of performing NLLSQ point mass fits to spherical harmonic representations of the Earth's gravity 
field [12]. Although these NLLSQ fits were, and are, of very high quality, the resulting NLLSQ problem itself 
is a very difficult one and it is perhaps best skipped whenever possible, unless the primary intent is to seek 
out a challenging problem (which is a worthy goal in-and-of itself). 

One recurring problem encountered in practice is that of implementing the analog of a regular gridding 
on the surface of a sphere. For example, severe angular distortion enters when a straight-forward evenly 
spaced angular grid is used, since the resulting sample points near the poles or z-axis are obviously closer 
than might normally be considered acceptable. While these unwelcome variations obviously do not occur 
for gridding on a plane, they are an integral part of the geometry of the sphere and cannot generally easily 
be overcome. The only way to completely overcome this problem is to place points at the vertices of some 
specified regular geometric solid figure (for example, at the vertexes of a cube), which fixes the number of 
interpolation points; however, this is not only very restrictive but more time consuming to implement for 
realistic general cases than it might at first appear. {Other strategies exist. One to consider might be that 
of placing a number of equal (repulsive) charges on a sphere and ascertaining the resulting configurations 
on the assumption that, in some sense or other, the resulting final equilibrium configuration will be that of 
a set of points that have the largest spacing possible. The theory of minimum energy point configurations 
on a sphere has its own extensive literature (see [6] and references therein), but even given this pre-existing 
literature there are greater implementation complications than one might want to deal with and, moreover, 
it is unclear whether this is really a general suitable point selection strategy to adopt for those not already 
familiar with the relevant literature in this area. A second strategy to consider might be one based on 
spherical i-designs — see for example, the first paragraph of [8\ and the citation given — but again it is unclear 
whether this is really a good strategy to adopt for someone who is not already familiar with the relevant 
literature.} For deep point source placements the issue of obtaining a more-or-less regular angular separation 
(i.e., one where the points are distributed somewhat uniformly over the entire surface of a sphere) is an 
important one. Specifically, from an examination of the values of the elements of the matrix T given by 
(gni), it is clear that as Pk' approaches Pk, Tkk' — > 1 (of course, Tkk' < 1 must hold for all Pk' ^ Pk)- This 
means that if two of the sample points on the sphere are close together, then the resulting C^^ will be large. 
By direct experimentation using irregular spacing on a sphere of radius Rp and comparing these results with 
results presented below, the reader can see for him or herself that somewhat regular spacing is very desirable. 
The question, then, is how to do this. An algorithm is introduced below in conjunction with Table [2] that 
makes a very good trade-off between ease of implementation and spacing regularity for the resulting grid. 
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Obviously the issue of point source spacing in the horizontal or tangential direction is the same for both the 
spherical exterior and interior problems, hence only the spherical interior setting will considered in the rest 
of this section. 

Having temporarily settled the gridding requirements for ip, next consider the issue of how to conveniently 
mathematically characterize the suitability of /, with regards to the overall angular separation of the spherical 
grid used. First, as alluded to above, consider the elementary well-known aspects of sampling along a straight 
line or of a 2-dimensional gridding in the plane. From the 1-dimensional theory of Fourier series and digital 
signal processing, the relevant concept for determining the appropriate spacing of an even grid of data is the 
concept of Nyquist frequency. Also the corresponding 1-dimensional criteria for smoothness of / is simply 
the highest allowed non-zero frequency term in a Fourier series expansion of /, which can be characterized 
by some integer n. The 2-dimensional Fourier series generalizations of this idea for a plane are obvious. 
For the surface of a sphere, the relevant concept might seem to be the maximum degree and order of a 
spherical harmonic expansion of / since a spherical harmonic expansion is the direct spherical analog of a 
2-dimensional Fourier series expansion for the plane. It is clear, however, that <j) and 9 are treated very 
differently in spherical harmonics since the undulations of say P„(cos 0) cos to0 correspond to very short 
distances at the poles for large n and m values; however, for small values of n (and thus m) this angular 
distortion is not quite so marked (i.e, for n and m less than about 7),. Thus for cases where n is small (or 
for cases where there is no dependence) , one can attempt to directly carry over the concept of a Nyquist 
frequency from lines in or planes in to spheres in R'^ and compute a Nyquist frequency based on 
cosmcj), for the maximum m, (or on cosnO for the maximum n) that may serve as a reliable guide for the 
sampling of / on a sphere [i.e., sampling of f{Rp, 9, (f>)]. It is clear that this spherical Nyquist criteria, which 
is characterized by the minimum great circle over which / varies appreciably, is the more meaningful one, 
rather than a direct characterization in terms of maximum allowed degree and order. For low degrees and 
orders the reader can easily test out the utility of this criteria for him or herself as follows. The idea behind 
this testing procedure is that it is a simple matter to distribute points evenly on an R^ unit disk and simple 
R^ harmonic test functions can be used to first see and then explicitly check out the role of the Nyquist 
frequency in this case; moreover, this simple testing algorithm need not be limited to R^ cases since any / 
that is a function of only x and y and is harmonic in R^ is also harmonic in M^. Looking slightly ahead, a 
more-or-less regularly spaced grid can be computed from the algorithm given in conjunction with Table [2] 
Suppose that such a semi-regular grid has been implemented, then consider the following simple harmonic 
functions in the interior of an unit disk: 

Him P^' sinM^ and H^m ■= P^^ cos AM (50) 

where p = \J + y"^ < 1, 'd — tan^^y/x and M is an integer. Consider the natural extension of these 
functions to the interior of a unit sphere, where the R^ x and y-axes are mapped into the R'^ x/y plane. In 
this case dHiM/dz — and dH2M/dz = so V^i?iM = V^iJ2M = for the interior of an R'^ unit sphere. 
Since an assumed grid of Pk values specifies the form of ^(X), it is a simple matter to fit to either Him 
or H2M for various integer values of M. Then, as the reader can then directly verify, the quality of such 
fits rapidly decline as the "spherical Nyquist sampling criteria" is violated. It might be argued that this is 
due to the special orientation that Hi m and H2 m have with respect to the 2-axis, but it is a simple matter 
to rotate the coordinate axes of either the sources or the iJjA/'s to show that this is not true. To further 
test the validly of this "spherical Nyquist sampling criteria" the reader can use a full spherical harmonic 
expansion to some appropriate degree and order in place of the Hj m^s if he or she desires, but the results 
should be obvious. 
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Specification of / and 



Since a discussion of most of the side issues is complete, consider the specification, in turn, of the actual 
forms used for / and ^ in testing. Fits to two different basic reference functions were performed: / = Fj for 
i — 1, 2, where 

• F2 := frPi(cos6l) - |r3p3(cos6l) + ^r'^ P^(co^d). 

Here the P„ are the usual Legendre polynomials [i.e., Pi (a;) = x, P3 = (5x^ — 3x)/2 and P5(a;) = (63a;^ — 
70x^ + 15x)/8]. [For consistency with earlier notation, where jo was used to denote an evaluation at the 
origin, and for the reader's convenience, the symbol Fj is used instead of /j, since fj is reserved to mean 

Both Fj 's were chosen because they are smooth functions and thus the assumed form for can reasonably 
be expected to work well with a small number of point sources. While the form of Pi was chosen mainly for 
its simplicity, the motivation for the choice of the form of P2 can also be easily understood. Thus let \Ja(^) 
denote the potential function that is harmonic for r < a and that is specified by the boundary condition 
/7a(a, 0) = 1 for 7r/2 > > and by /7a{a, 0, ^) = —1 for tt > > 7r/2. P2 then corresponds to the first 
three terms of the spherical harmonic expansion of Ua{X), with a = 1 [11]. While P2 is smooth at r = 1, 
it is does have a somewhat abrupt change in transitioning from 9 < 7r/2 to 9 > 7r/2. Notice that when 
a — 1 the full expansion for Ua is not harmonic on the surface of the unit sphere due to the discontinuity at 
9 — 7r/2 and thus it is not an allowed function (since > 1, only functions that are somewhat smooth on 
the unit sphere itself can be accurately modeled) . Furthermore, partials with respect to z are not continuous 
at z = and x"^ + y'^ = 1. 

Next, consider the actual specification of (p. The scalar point source basis function sets can be fixed by 
specifying the location of interpolation points (Pfc) and whether a constant basis function term is present or 
not in the fitting form. First consider only implementations without a constant term. Further, as discussed 
above, suppose that the interpolation points form a grid of points on some sphere of radius Rp < 1. Then, as 
previously noted, from the collocation property one might reasonably expect that / itself will be accurately 
matched over the entire sphere of radius Rp < 1 from the uniqueness of Dirichlet boundary conditions; 
furthermore, one might expect that / and ip will match reasonably well over the part of the interior region 
specified by Rp > r, due to the matching properties of the II norm itself. Thus, it is a reasonable hypothesis, 
pending numerical testing, to assume that ip and / will also match over all of the interior region, with the 
largest excursions occuring near the surface of the unit sphere. On this theory the basic configuration tested 
will thus be that of a set of interpolation points positioned on a sphere of radius Rp < 1. This is done not 
only for convenience, but also for conceptual clarity since there is otherwise no obvious strategy, at present, 
for efficiently placing point sources so as to produce matches to /'s that share common characteristics. Also 
as noted above, while it is clearly desirable that this set of interpolation points be uniformly spaced over this 
sphere, this cannot be done in general exactly for an arbritrary number of points; moreover, it is a nontrivial 
matter to implement the particular cases where such a regular spacing is possible. To address this problem 
a very simple algorithm is given in the next paragraph that produces a sufficiently regular spacing of points. 
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Symbol 


Input Ng 


Nk 


Mill. Ang. Sep. 


Max. Mill. Aug. Sep. 




8 


58 


25.71 Dcg. 


27.96 Dcg 


CP2 


36 


1542 


5.14 Deg. 


5.30 Deg 



Table 2: Point Configurations, Cp.{Rp), for Sphere of Radius Rp 



Specification of Interpolation Point Gridding 



The interpolation (and source) placement algorithm places points on various rings and is based on main- 
taining a minimum angular separation of points while simply adding as many points on each ring as can be 
added without violating this spacing constraint. Since these points on the sphere of radius Rp are deter- 
mined by a grid of 9 and ((> values, which are fixed, in turn, by an array of delta spacing values A9 and A0, 
it is convenient to define this grid of points on a unit sphere and then simply rescale them to the appropriate 
value of Rp, as needed. 

First, the angular separation, Pkk', of two points Pk and Pk' will prove to be a very useful concept and 
it is defined by 

Second, the algorithm to fix this grid of 0's and 0's is uniquely determined by the input number Ng, which 
specifies the number of rings parallel to the x/y-plane and whose z coordinate locations are fixed by the 
different values of 9 (i.e., z ~ cos 9). Let the indexes n and m be associated with the various 9 and <j) values 
respectively. Then A9 := tt/^Nq — 1). The first "ring" (n = 1) consists of the single point (0, 0, 1)"^ along 
the positive z-axis and the last ring [n = Ng) consists of the single point (0, 0, —1)^ along the negative 
z-axis. For 1 < n < Ng the ring locations along the z-axis are given by dn ^ [n — 1)A6'. For each value of 
n such that 1 < n < Ng, there are associated (j) values that fix the points on ring n. For a given n, let A(/)„ 
denote the angular spacing for ring n. For a given A(/)„ the values of cf) along the ring will then be given by 
4>nm '■= {m— 1)A0„, where m = 1, 2, 3, . . . , N^ and A0„ :— 2it/NcI,. This leaves only N^ to be determined, 
which is the maximum integer that will ensure that all points on the ring obey the minim angular separation 
requirements. To compute N^ simply test out successive values until the inequality Pkk' < A0 is violated 
for two typical adjacent points on the ring and then set N^(n) to the last value for which this inequality 
held. Processing the various rings in sequence and then rescaling to the appropriate Rp, this algorithm then 
completely specifies the Pk by (r, 9, (p) — {Rp, 9n, 4>nm)- 

For Ng — 8, this algorithm gives Nk — 58. This configuration of points is labeled Cp^{Rp) and is 
completely specified by an input value of Rp. For simplicity this configuration is the one that is used to 
define ip. For Ng = 36, this algorithm gives 1, 542 points. This configuration of points is the one used for 
evaluating the quality of the resulting fits and is labeled Cp^ [Re)- It is completely specified by an input for 
the evaluation radius. Re- For any of these point configurations there is a minimum angular separation of 
distinct points and a maximum minimum angular separation of distinct points, which together give some 
indications of the uniformity of the distribution of points on the sphere. The minimum angular separation 
is defined by the minium of [3k k' for all k and k' ^ k. For each k, let jk ■= minimum of f3kk' for all k' ^ k. 
Then the maximum minimum angular separation is given by the maximum of jk for all k. These angular 
separations for Cp^{Rp) and Cp2{Rp) are given in Tabled 
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Case 
Number 


/ 

Choice 


Rp 

ipk) 


if at 
Origin 


Fit Errors at = 1 


Fit Errors at = 1/2 


Standard 
Deviation 


M cLxi mum 
Magnitude 


St&.nda,rd. 
Deviation 


Maximum 
Magnitude 


1 




/zt 




.luUX lU 


.0141: X lU 


970 V 1 c\~^ 
.Z 1 U X ±u 


.^yo X lu 


9 


^1 


i/Z 


1 nn V 1 n~6 


'^fi9y 1 n~3 

.0\J^ A i_U 


. XUO A J.U 


. XOZ X xu 




•3 
O 




1 /I n 


9zLzl V 1 n^i''' 


QS'^ X 1 n-s 

. C/OO A 1 U 


407 V 1 0"'' 


.O ( 4: A 1 U 




4 


-f^ 1 


1 

100 


n Dx 1 n° 


.947x10-" 


.387x10-13 


.370x10-1° 


.151x10-1^ 


5 


i^l 


9 

1000 


.197x10-29 


.511x10-1^^ 


.208x10-13 


.200x10-1° 


.826x10-1° 


6 


F2 


3/4 


-.389x10-3" 


.519x10-1 


.131x10° 


.827x10-3 


.127x10-2 


7 


F2 


1/2 


-.277x10-31 


.931x10-2 


.371x10-1 


.304x10-"^ 


.587x10-4 


8 


F2 


1/10 


.473x10-29 


.135x10-4 


.593x10-'' 


.265x10-'^ 


.114x10-° 


9 


F2 


1 

100 


0.0x10° 


.137x10-8 


.589x10-8 


.277x10-11 


.114x10-1° 


10 


F2 


9 

1000 


.290x10-23 


.109x10-8 


.427x10-8 


.311x10-11 


.107x10-1° 



Table 3: Fit Results to f ~ Fi or F2 with 58 Interpolation Points 



Numerical Results and Discussion 



Interpolation results for fits to Fi and F2 using Cp^ at different Rp's for the configuration of ip are 
presented in Table [31 The first three columns of Table [3] are self-explanatory. The fourth column gives the 
value of ip{6) = (fiiO) - Fi(0) = .^(0) - ^2(0), since Fi(0) = and ^2(0) = 0. The fifth and sixth columns 
then give statistical parameters for the fit based on Cp^ evaluated at Re — 1- Like results for Re = 1/2 are 
then given in columns seven and eight. 

Additional information about these fits is given in Table |4l where for convenience the first three columns 
are repeated. In Table |4] the fourth column is the collocation check number, which is the maximum of 
\f{Pk) ~ ^^Pk)\ for all fc, and gives some measure of numerical roundoff. The fifth column is the condition 
number as previously defined. The sixth and seventh columns gives the statistical parameters for the source 
strength (gfe), while the eighth and ninth columns give the fit error parameters evaluated ciX Re — 1/4. 

There are several significant observations that can be immediately made from Tables [3] and [U which tie 
into previous discussions: 

(a) Placing the interpolation points on a sphere of radius [Rp < 1) rather than distributing them through- 
out the whole interior of the unit sphere seems to work well for the tested functions. 

(b) Both F\ and F2 are very smooth functions. Small choices oi Rp, large C^'s and accurate fits occur 
together for these test functions. As previously discussed, while small values of Rp and large C#'s 
are unavoidably linked for smooth functions, for irregular functions better results generally occur for 
larger values of Rp and more interpolation points. 

(c) From the fourth column of Table 2] it is apparent that all of the fits satisfy the collocation condition, 
but as one might expect a small degradation occurs for large condition numbers due to the onset of 
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Case 

# 


/ 

Fit 




Coll. 
Check 


C-u- 


Source (g^) Statistics 


Fit Errors at i?^ = 1/4 


Mean 


Staiidcird 
Deviation 


Sta,iida.rd 
Deviation 


Maximum 
Magnitude 


± 




J/ 4 


1 R V 1 n~32 


.UU A lU 


. 1 ( A lU 


01 1 V 1 
.y ± ± X lu 


.yoo X iu 


.zuy X lu 


9 






1 '^v 1 n~32 


.4:?? A ±U 


i^i V 1 n~8 


'^1 n V 1 0^ 

.O XU A iU 


. J. X L\J 


.ZrO<J X J.U 


o 




1 /I n 


71 V 1 n~3i 

. / 1 A lU 


so V 1 f)!^ 
.OU A lU 


49 V 1 n~18 


. OOO A ±u 




.4:Uiy A ± U 


4 


Pi 


1 

100 


RSv 1 0~2^ 


SDx 1 0^3 


.98x10-3" 


.388x10^ 


.145x10-1^ 


.592x10-1® 


5 




9 

1000 


.14x10-28 


.43x1034 


.28x10-29 


.532x10^ 


.781x10-19 


.333x10-18 


6 


F2 


3/4 


.18x10^32 


.66x10^ 


-.90x10-32 


.253x10'^ 


.107x10-3 


.270x10-3 


7 


^2 


1/2 


.43x10^31 


.49x10^ 


-.10x10-32 


.267x10^ 


.170x10-^ 


.330x10-5 


8 


^2 


1/10 


.28x10-28 


.80x10" 


.87x10-30 


.402 xlO^ 


.519x10-1° 


.204x10-9 


9 


^2 


1 

100 


.15x10-23 


.80x10^3 


.45x10-24 


.402x10" 


.614x10-14 


.219x10-13 


10 


^2 


9 

1000 


.10x10-22 


.43x10^4 


.37x10-23 


.756x10" 


.120x10-13 


.404x10-13 



Table 4: Ill-conditioning Indicators for Table [3] Fits 



numerical roundoff. 

(d) C^'s are sufficiently large in Table [3] that the issue of finding ways to naturally reduce them arises. 
One way will be considered below, but other ways will not be addressed until subsequent articles. 

(e) As noted earlier, and as indicated in Table [31 large condition numbers are associated with large excur- 
sion in the source parameters. 

(f) Since the results of Table [3] indicate that the approach works well for smooth functions, this suggests the 
strategy of residual fitting: i.e., first fitting the smooth part of a harmonic function and then removing 
this part and then performing a fit the remaining part, etc. The final result can be reconstructed as a 
sum of these two (or more) parts. (As just noted, these and other similar strategies will be addressed 
in another venue.) 

Points related to the above items dealing with large C^'s, as well as to the ideas discussed earlier, can be 
found in, for example, the pair of articles [TB] and [3]- 

To bring other points to light, additional testing was done based on the results of Tables [3] and 2] (most 
of these tests are modifications to the basic configuration of test Case 3). First, to see how an additional 
sphere of point sources would effect the results, the configuration Cpj (1/2) U Cpj (1/10) was run and, as one 
might expect, resulted in a slight improvement over the results of Case 3: A standard deviation of .202 x 10-^ 
resulted, but with a larger associated of .43x1022. 

Second, there are three ways to account for the occurrence of constant terms in /: 

(A) Use an unmodified basis function configuration to sec if constant terms are naturally handled properly. 

(B) Add a constant basis function directly to ip. 

(C) Just subtract off the constant term from / and then add it back on. 
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(D) Add a point source term near the origin. 

To test out the how fits work when a constant bias is added to /, the following fits were done using approaches 
(A) and (B), which both seemed to work well: 

• Case 3 results were replicated with Fi replaced by F3 := Fi + 1. The quality of the fit was unchanged, 
which indicates that constant terms can be naturally handled if the point sources are reasonably 
positioned. 

• The Case 3 fit to was rerun with a constant basis function added to ip. An improvement in the 
quality of the fit resulted (with a standard deviation of .304x10^^ at Re = 1) since the added basis 
function implies more fitting degrees of freedom; however, the was somewhat larger (.11x10^^). 

Approaches (C) and (D) were not tested. 

Third, to test out the effects of normalization of basis functions, the following experiments were run: 

• The results of Tables [3] and |4] were rerun using normalized basis functions and the condition numbers 
were unchanged. Furthermore, aside from expected numerical changes (in, for example, columns labeled 
"ip at the origin", "Coll. Check" and the Mean value of q^) the results were almost identical. 

• The configuration Cpj (1/2) U Cpj (1/10) mentioned above was rerun with normalized basis functions 
and had an improved C# (of .29x10^^). 

• A fit using the ip specified by Cp^ (1/10), with a constant term added yielded a larger C# (of .15x 10^'^). 
This indicates that normalization does not necessarily always improve the C^^^; however, as previously 
noted it is definitely called for when mixed dipole and point source fits are performed due to the 
disparate nature of the physical scaling involved. 

A final word of caution is in order here. Although the quality of the match to the sampled values of a 
smooth function may be very good in the presence of large C^f^'s, small deviations due to either a rougher 
function than expected or to the existence of measurement sampling errors can be of special concern. In 
particular, it is obvious that the propagation of errors always deserves special attention when large C^'s 
are present. As the reader can verify for him or herself directly by perturbing individual components of the 
A vector, any miss-match here, even a very small one, can result in an extremely poor quality fit and, in 
fact, a small error at some interpolation point(s) may well result in a fit error at certain field points X that 
is much larger than the value of |/(A")| itself, even though the fit still matches at all of the interpolation 
points: That is, in the presence of large C#'s, if measurement errors are present then |</5(X)| >> |/(^)| for 
some values X ^ Pk may well occur. One obvious strategy is to increase the number of interpolation points 
and place them closer to the surface of the sphere. There are various other strategies and implementation 
procedures that can be considered, but as mentioned in Section [1] errors and ways to treat them are not 
considered in the present article. Finally, as discussed above, if / is unsuitable or if the sampling points are 
inappropriately located |'y3(A')| >> |/(A')| for some values X ^ may also happen due to other causes 
and this is one of the main reasons that application of the text approach have been limited here to modeling 
problems where these possibilities can be explicitly checked for by means of direct comparisons. 

6 DIDACKS MFS/BVP Connections 

At the end of Section [1] it was pointed out that the technique presented here can be profitably used in 
conjunction with MFS BVP approaches for certain problems, thus the use of the current method in conjuction 
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with BVPs, and especially the MFS approach, is singled out for further discussion in this section. Only a 
preliminary analysis of DIDACKS connections to MFS approaches will be undertaken here. 

Prior to considering MFS approaches themselves, consider how most of the numerical test cases given in 
Section O can be reinterpreted as solutions to spherical BVPs for some interior region of 51 [qj . First recall that 
the general source basis function placement strategy used in Section [5] was to place sources on an exterior 
sphere of radius l/Rp > 1, which corresponds by (jlOp to placing interpolation points on a sphere of radius 
Rp < 1. From Theorem [l] since the fitting function ip matches the function to be fit / at the interpolation 
points, when / satisfies an appropriate spherical Nyquist criteria, as discussed in Section [Sj (p and / match 
over the sphere of radius Rp and (p is a. solution to the Dirichlet BVP for r < Rp. Case 2 and 7 of Table [3] 
correspond to solutions of Dirichlet BVPs for the interior of a sphere of radius Rp ~ 1/2, where if matches 
Fi over r = Rp for Case 2 and ip matches F2 over r = Rp for Case 7. The last two columns of Table [3] give 
the corresponding mismatch over the boundary of this sphere of radius r = Rp, while the last two columns 
of Tabled] give the corresponding mismatch over an interior sphere of radius r — Rp/2. From the discussions 
in Section 3] and Appendix A these BVP solutions also satisfy the following two norm minimization criteria: 

(A) 11/ — (^lljj is minimized. 

(B) ||<^||jj is minimized. 

While these norm minimization conditions hold for the whole of f2[o] and not just that part inside the sphere 
of radius Rp, they still are still clearly very significant properties of the BVP solution (p. While it is possible 
to analyze the various mathematical implications of (A) and (B) and then derive related inequalities for 
the subregion r < Rp; for our purposes here a direct analysis of (A) and (B) will have to suffice. Thus 
notice, for example, that statement (B) makes a fairly strong assertion about the harmonic extension of ip 
from the region r < Rp to the whole of region il[o]- (For a- further discussion of the properties of harmonic 
extrapolation see Appendix B.) Finally, it is obvious that if ip matches / over the whole of fi^o] then it must 
match over the subregion r < Rp also: In fact, from the properties of harmonic functions, since the largest 
excursion in f — ip must occur at r = 1 one would naturally expect a much better fit over the subregion 
r < Rp than over the subregion Rp < r < 1, as seen in Tables [3] and [4] for the two cases noted above. Third 
and finally, as indicated in Table [4] the C^'s for most of the cases dealt with here are much larger than those 
normally encountered in MFS applications where, as noted below, many more sample points are commonly 
used. 

Several generalizations of the basic BVP approach just outlined are obviously possible. First, as men- 
tioned before, it is a trivial matter to rescale the results in order to handle interior problems for spheres with 
an arbitrarily specified radius. Second, since the various partials of / correspond to higher multipole fits it 
is easy to consider BVPs for interior spheres that have other specified boundary information than potential. 
Third, as previously mentioned, corresponding approaches to BVPs over unbounded domains can be based 
on the I norm. Fourth, the above ideas, including the discussion of norm criteria, apply to other DIDACKS 
settings, such as to the interior and exterior of units disks in the real and complex plane [20] . Fifth, since 
DIDACKS theory does not otherwise constrain interpolation point locations, it is quite possible to consider 
a like approach to BVP problems with more general shaped boundaries. This last topic will be taken up 
below. 

Next specifically consider how the MFS approach relates to the basic BVP technique just outlined. In 
what follows, it is assumed that the reader has at least some familiarity with MFS approaches — the reader 
who is not familiar with MFS approaches can consult any of the outlines that are normally contained in the 
front of current articles dealing with MFS techniques, such as [7], [21], [22], [18] or [15j (all these references 
also give additional references). Bogomolny [3] is a frequently referenced article that contains probably the 
first mention of the confluence of good solutions and large C^'s for MFS-like approaches. One normal way 
that MFS is applied to the interior or exterior of a sphere where Laplace's equation holds and Dirichlet 
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boundary conditions are specified at some radius, say r — Rp, is by simply postulating a linear combination 
of the form Lp := X^fc^i Qk/^k and then determining the associated source strengths qfe S M from a forced 
collocation point matching condition analogous to (j38p . which is a strategy called direct collocation by 
Saavedra and Power [52] and that is discussed at length in Appendix A in related contexts: 

^{Pk) = f{Pk) for fc = 1, 2, 3, . . . , iVfe . (52) 

Here the interpolation points are on the surface r = Rp so Pk = Rp and the source points are located 
on some other spherical surface, say r — Rx' , so that r'^ := ~ Rx' ■ For MFS problems with 

this geometry, it is usually further assumed that Pk/Pk — ^k/^'k^ which means that for any Rp and Rx' 
when the problem is rescaled in an appropriate way (jlOp can be assumed to hold. When this geometry 
holds then the fit corresponding to is a DIDACKS scalar point source fit that, as discussed above, is 
a solution to the associated Laplacian EVP. For this geometry, given this direct equivalence of MFS and 
DIDACKS approaches it follows immediately that the minimum norm conditions (A) and (B) outlined above 
apply and that they have the indicated interpretation. Furthermore, other conclusions can be drawn from 
discussions of DIDACKS theory given previously, such as the explanation for the confluence of large C#'s 
and well behaved solutions. Finally, even after taking the implied correspondence between DIDACKS fits 
and MFS fits into account, observe that in the MFS literature there are several differences between the way 
fundamental solution basis function fits are commonly implemented and the way they were implemented in 
Section [5l First, many more fundamental solution basis function are generally used than indicated by the 
configuration Cp^{Rp) of Table [2l Second, the source locations indicated in Table [3] are generally much 
deeper than is normal in the MFS literature. For example, for spherical interior problems cites a general 
rule-of-thumb that the magnitude of the source point vectors should be five times the magnitude of the 
interpolation point vectors or = 5Pfc. From (fTO|) this means that = 5 or that Rp ~ \/5 ~ .45, which 
corresponds closest to Cases 2 and 7 in Table [31 

In the literature MFS approaches to C or unit disks are also common. For complex or real unit disk 
geometries the same conclusions can be drawn, but complications may enter since sometimes the interpolation 
points are offset by a certain angular amount from the direction implied by the corresponding source vector 
so that the relationship for the plane corresponding to ([TOj) fails to hold. Besides the simple direct collocation 
approach outlined by (|52p and the complication just mentioned, there are two other common ways that MFS 
techniques are applied. In thee cases, it is still often possible to draw significant implications from DIDACKS 
theory about the associated MFS solutions; however, these discussions are outside the scope of the present 
article. There is one area that is of more immediate concern — namely the question of how general the 
boundary surface can be with regards to DIDACKS approaches. Notice that even when the interpolation 
points are placed on a different shaped interior surface in Sl[o], norm conditions (A) and (B) still hold and 
while they still have considerable general significance their physical meaning becomes less clear. 

The possibility of applying DIDACKS theory to general bounded domains is easy to understand. Consider 
a (possibly dense) set of interpolation points specified on some closed star-shaped surface S (c.f, Definition 
B.l), which lies solely within some spherical domain. For a suitable given spherical domain, DIDACKS 
theory then prescribes the locations and strengths of a set of point sources that will replicate the specified 
information on E. The sources are located on some corresponding surface E', which is assumed to be wholly 
in the exterior of the given sphere. Clearly if the interpolation points are properly chosen then the DIDACKS 
fit yields a function that is harmonic inside the region bounded by S and that matches the boundary value 
information specified on it. Since the underlying fitting process reconstructs the BVP as a Poissonian inverse 
source problem and these particular inverse source problems are notoriously ill-conditioned, some restrictions 
on the shape of E and on the locations of the interpolations points, which are on S, are clearly called for. 
For example, for a given choice of sphere due to the fact that the source locations themselves are directly 
specified by this BVP information, the resulting source locations may be entirely inappropriate. This clearly 
opens up a number of issues related to both the applications and theory of 
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Additional issues also merit further examination here. For purposes of discussion, consider only the 
exterior problem for unbounded regions, which have a single bounding interior surface. First, even for 
BVPs with fairly irregular geometries there may be no real problem if the underlying field function to 
be matched displays very benign behavior and the source locations are appropriately chosen — that is, if 
various interpolation points on the boundary surface are naturally related to the given (benign) attenuation 
properties of the harmonic function that is to be reconstructed. (It is worth emphasizing that even in this 
situation the choice of implicit sphere used for positioning the sources must be made fairly carefully.) For 
optimal solutions, some problems require both fairly irregular shallow point source placement as well as 
deep source placement. For such problems, while it may be that the natural interpolation points fall on an 
appropriate prescribed BVP surface that has suitable boundary conditions, this clearly would be very much 
the exception rather than the rule. As indicated in Appendix A, when the properties of / are unsuitable 
and the sources are badly placed the resulting error excursions may even be many times larger than the 
maximum of |/| itself in much of fl. 

Second, for nearly spherical geometries the problems encountered may be minimal, so fairly direct schemes 
can be used to overcome any residual problems. In geophysics, for example, there are fairly standard 
techniques for addressing gravitational problems that can be adapted to some other scenarios [10] when 
the geometries are not to irregular. The techniques used there include always subtracting off a standard 
ellipsoidal reference model from the total field to obtain what is called gravity anomaly (or a related quantity 
called gravity disturbance). Often a further subtraction, which is called remove-and-restore, is carried out 
that results in a localized field. At several other places in this article the advantages of such residual fitting 
techniques for DIDACKS applications have been pointed out. General issues associated with sources near 
or slightly outside the harmonic region of interest are also common to physical geodesy; for example, over 
land where mountains protrude above the reference ellipsoid, reduction techniques are employed to obtain 
equivalent sources that do not extend above the ellipsoid — for additional points see [TUl HT] . 

Third, and finally, there are certain geometric consequences of the DIDACKS relationship (fTU]) between 
field and source points that must be considered. For the sake of argument, consider an exterior harmonic 
BVP problem for some surface S that has been positioned and scaled so that a unit sphere that is constrained 
to be inside S is a reasonable candidate choice for the source region. Suppose that the sample point 
P e S C M'^ is located at some peak of S, then from (fTO|l the corresponding source placement will not be 
closer to the surface as one might expect, but further away and, in fact, deeper than the other source points. 
This geometric effect is labeled the reverse-dimple effect. For nominal DIDACKS fits this reverse-dimple 
effect is quite natural since sample points that are further away will be markedly attenuated and thus deep 
sources are much more appropriate; however, for BVPs, or more specifically BEM or other source coating 
techniques, one would nominally expect to find sources at and around excursion points like P. Hence, the 
reverse-dimple effect poses a conceptual problem in adapting DIDACKS theory to these BVP approaches 
for non-spherical geometries. While there may be various ways of overcoming these serious conceptual issues 
for non-spherical problems, any real progress in this area is outside the scope of the present work and must 
await the efforts of others. (Keep in mind, however, in the end that the underlying DIDACKS formalism 
makes no restrictions on the locations of sample or source points, aside from the obvious requirement that 
the be in the correct region. Also, in passing, notice that in the appendices P will not take on the special 
meaning that has been assigned to it in this section.) 
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Appendix A 



The Collocation Minimum Norm Property 

In Section IH the collocation minimum norm property for scalar point sources over r2[o] equipped with the 
II norm was stated by ([38)) and in a slightly more general way by Theorem [1] This appendix entertains the 
possibilities of various other generalizations. It also discusses the specific ramifications for RKHS based fits. 
Much of the exposition here has a work-in-progress feel to it and the bulk of this appendix is not aimed at 
the general reader, but rather to those either engaged in or contemplating research in the areas of RKHS 
theory, Bergman kernel theory, general kernel theory or in the areas dealing with the interaction of general 
collocation kernel theory and operator theory. 

For convenience this appendix is largely independent of Section [H so a restatement of ([38]) for a slightly 
generalized setting is given first. Thus, here the minimum norm property means that of all admissible 
functions, g such that 

g{Pk) = Cfc for fc = 1, 2, 3, . . . , iVfe , (A-1) 

where the Ck are arbritrary constants, a DIDACKS scalar point source fit with coefficients qk specified by 
(l2|) or a proper RKHS fit is the function of smallest norm: \\g\\ > \\ip\\. The general plan of approach will 
be to explore the limits of what can be reasonably obtained in the way of possible generalizations, before 
giving formal results. Thus (|A-ip will be analyzed and restated in slightly different ways as the section 
progresses. Although the resulting exposition may be somewhat circuitous, various subtleties will emerge as 
the section unfolds and the reader should finally end up with a better overall understanding of the relevant 
issues involved than might have resulted from a more direct exposition. 

The corresponding minimum covariance norm result for various norms over VLq is a cornerstone of GC 
|17j . which as previously noted, is often called linear least squares collocation. Related results are known 
to hold for other reproducing kernel contexts, such as analytic function theory and functional analysis. For 
a general RKHS setting this minimum norm condition can be stated as: Of all those functions that match 
some specified function (/) at Nk points, a linear combination of (symmetric) reproducing kernels is 
the one that has the smallest norm. This will be labeled the RKHS minimum norm property (this section 
presupposes a rudimentary understanding of RKHS theory — see, for example, [Tl[16l|23].) There are several 
relevant issues that need to be considered in association with the conditions required for this RKHS minimum 
norm property to hold, which are not explicitly addressed in standard references. Furthermore, not only is 
there no existing reference known to the author that discusses these subtleties, but many references seem to 
be, at best, slightly misleading. Some additional attention and special care is needed here since this result 
must be extended to DIDACKS theory. 

This appendix is fairly long since a somewhat thorough, but informal, analysis is given. Prior to Theo- 
rem [3] most of the discussion is an effort to establish a proper framework and to determine the appropriate 
limits of its generality. In the end, the easily codified results are directly linked to the definitions of properties 
of the class of admissible functions and of the class of operators that are allowed to act on them. With re- 
gards to relevant additional trait that any allowed family of admissible functions is assumed to satisfy, while 
other characterizations of it are considered prior to introducing the idea of pointwise independence, this is 
the trait that is deemed most appropriate in the end due to the various considerations brought up. With 
regards to the relevant characterization of linear differential operators, due to other types of considerations, 
a slight generalization of the one considered in Section |4] is finally settled on and introduced by Definition 
A. 4. These two concepts immediately yield Theorem [71 although it is, in reality, only a slight generalization 
of the results of Section [4) On first reading the reader may wish to skip directly to the Formal Result given 
at the end of the appendix and in so doing the exposition gains a certain measure of precision, cohesion 
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and closure, but the reader will do this at the expense of a deeper overall understanding and knowledge 
base. The intervening material here before this last part of the appendix is not especially well polished or 
organized and includes some topics, such as geographic factors, that are entertained primarily to simply 
gain a different perspective on existing material. Finally, what is labeled the collocation minimum norm 
property or simply the minimum norm property might just as well have been labeled the RKHS minimum 
norm property except for the fact that a more generalized setting is addressed here. 

Preliminary Comments 



As previously noted in Section [H there are essentially two ways to approach issues associated with 
minimum norm properties: Type I proofs, which are based primarily on various properties of the fit and 
its norm \\{p\\. Type II proofs, which are based primarily on the properties of the function to be fit, /, as 
well as its associated function space. The former is the one taken by proponents of GC, while the second 
is the one frequently taken by mathematicians. The informal analysis undertaken here starts by examining 
Type I approaches in detail before finally considering Type II approaches, which in turn leads back to 
considerations associated with Type I approaches. As noted below, there are issues with both types of 
approaches. Although Type II approaches seem to be safer, one historical problem with Type II approaches 
is that results are generally associated with a functional analysis setting, where it is obvious that each point 
of / can be varied in dependently and this attribute is taken for granted. This may not be true even when 
very natural restrictions are implied by the nature of the family of admissible functions. For both types of 
approaches there seems to be a need for a commonly unstated requirement that Nk < Nf, where Nf is the 
dimension of the span of the family of admissible functions. Otherwise, for example, if Nk = Nf a trivial and 
uninteresting situation results where there is at most only one function that satisfies (jA-ip . so that (fi = f 
and thus ||<y3 — /|| =0 and \\ip\\ — \\f\\. 

For Type I approaches various other issues are involved. Hille seems to be fully aware of all the relevant 
issues involved and gives a clear and rigorous proof of this RKHS minimum norm property for a somewhat 
limited case (i.e., for N^ — 1 only in the one dimensional complex setting) ^ p. 333]; however, as noted 
above, there appears to be no reference that discusses the relevant subtleties involved and simultaneously 
gives a clear, concise and rigorous discussion for general cases. The authors that are mindful of the issues 
involved generally seem to give either somewhat limited or specialized piecemeal discussions. (Furthermore, 
even for the restricted RKHS setting, a discussion of the collocation minimum norm property often does 
not appear in the mainstream literature at places where one might expect to find mention of it, such as 
[1] and [23].) Most, if not all, of the issues involved concern either the class of admissible functions or the 
invertibility of the associated Gram matrix [T in ^] or both. For RKHS theory the usual form of the 
underlying Gram matrix can be inferred from (jA-7[l . Let G denote this matrix in general. Since a Gram 
matrix is the Hermitian matrix of inner products of basis functions its, eigenvalues are greater than or equal 
to zero and either |G| > or |G| = 0. In keeping with current conventions, when |G| > holds G will 
said to be positive definite and if |G| = cannot be ruled out, then it will be called positive semi-definite. 
It is thus the case |G| = that causes problems, since it implies that the associated basis functions are 
not linearly independent. Historically, many authors, such as Aronszajn [l], used the term positive definite 
for what is labeled positive semi-definite here and presumably this variation in nomenclature has given 
rise to some of the confusion that has occasionally arisen in the past. As is often the case with matters 
pertaining to linear independence, there are greater subtleties involved here than one might at first expect. 
After the main results of this section are derived, a discussion of these matters is included in an attempt to 
provide the reader with some of the tools necessary to sort out the main issues involved on a case-by-case 
basis — especially with regards to the RKHS context. Hopefully, it should be apparent from the sequel that 
simply adding the explicit assumption of invertibility of the underlying linear systems to existing RKHS 
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minimum norm demonstrations is all that is needed in order to take care of all, or at least most, of the issues 
involved. With regards to RKHS collocation, although the discussion here may seem to be overly elaborate, 
the relevant points involved are considered well worth making. 



Point Collocation Without Differential Operators 



Before considering the complications that result from allowing point collocation or replication conditions 
that have associated operators, it is useful to consider the implications of a direct point matching condition 
like (jA-l[) . At the end of this discussion a simple way of describing most of the relevant points will emerge 
from an examination of standard RKHS fits. 

In what follows a general M" setting is assumed and in order to set this setting apart from that of the 
rest of the paper a new notation is used for it, but only where it is appropriate. Thus, at the risk of being 
repetitious, let Q now denote some appropriate region in M", z a point in Q and f{z) a typical member of 
some admissible class of functions defined over SI. It is assumed that H. is connected and does not consist of 
the entire space R" (so that if the class of admissible functions is harmonic, then f{z) will not necessarily 
be zero). In what follows all the basis functions and kernels will be assumed to be bounded members of 
the associated function space. Further for some relevant inner product and positive definite norm structure 
(11/11 > for all / 7^ 0), let the cost function to be minimized be 

^■.= \\f-^r = \\fr-2{f,^) + M\ (A-2) 

where, as usual, / is the function to be approximated and ip is the approximating function. Here ip is assumed 
to be a linear combination of fixed basis functions of the form 

ipiz)=J2t^kBkiz) (A-3) 

and these basis functions themselves are assumed to also be in the class of admissible functions. Since the 
class of admissible functions is assumed to admit linear superposition, a finite linear span of admissible 
functions yields an admissible function and so ip is an admissible function. For convenience, it is further 
assumed that f{z), Bk{z) and /Zfe are real valued. 

As before, the values of /ifc that minimize <& ~ $(Aife) are: 

Y,iBk',Bk)fik^{Bk', f) . (A-4) 

k=l 

which can be written more compactly as 



T^ = A (A-5) 



where T again denotes the matrix whose elements are Tk'. k = {Bk' , Bk), A. the vector whose element are Ak = 
{Bk, f) and fi the vector whose elements are /Xfc. Functional minimization processes with an underpinning 
inner product structure will labeled generalized least squares (GLLSQ) processes, so as to distinguish them 
from the more common processes where a sampled sum of squared residual errors is minimized. This 
standard sampling least-squares process will be labeled linear least squares (LLSQ). For a LLSQ fit using a 
combination of Nj^ basis function, if there is a sum of Nji squared residuals then in general Nj^ >> N^- 

Before proceeding to a consideration of the main theorem dealing with minimum norm properties, the 
following definition is useful: 
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Definition A.l The closed form generalized linear least squares (CFGLLSQ) condition is said to hold for 
an inner product structure and approximating form specified by (|A-3p when it is possible to evaluate {Bk' , /) 
in closed form and Equation (jA-Sp is solvable (i.e., exists). Further, the vector /i is assumed to be fixed 
by (jA-Sp . When the kernel form Bk' is a DACK and produces closed form expressions that are a function 
of a single vector argument for (Bj^,, /) [such as (jl7p ]. then it is called a point replication kernel. Also as 
noted above, / and Bk are assumed to be members of some class of bounded admissible functions so that 
11/11 < oo and \\Bk\\ < oo. 

Observe that Definition A.l only requires that each Bk be associated with some argument in . For 
convenience let tk € denote this argument. In what follows it is assumed that there is a one-to-one 
correspondence between the {Bk} and {tk}'- tk = tk' if and only if k — k'. Further, in the remainder of 
this section it is assumed that the values of /ife are determined by (|A-5p so that ip is a fixed function (c.f., 
Definition A.l) and that |^fe| < cx) so that \\ip\\ < oo. 

Definition A.l is quite general and it is worth considering what types of kernels are allowed. First, all 
of the DIDACKS examples considered in Section [2] satisfy the CFGLLSQ condition, including dipoles and 
higher order multipoles, but this fact will not be shown here. Mixed types of fits are also allowed since each 
Bk can be a different type of basis function. Second, reproducing kernels, K(x, tk') for tk' G f^, with an 
associated RKHS and norm are allowed if the associated exists, since by definition {K{tk' ,z),f) — f{tk'); 
moreover, in the real setting when a linear combination of reproducing kernels of the form 

^{z) = Y, l^kKitk, z) (A-6) 
fe=i 

can be fit to /, then (jA-4p can be rewritten as 

Nk 

Y,Kitk',tk)tik^ f{tk') (A-7) 
fc=i 

since Tk' k = {K{tk',z), K{tk,z)) = K(tk', tk) and K{tk', /) — f{tk')- Third, Dirac delta function forms, 
such as 6{x, tk') = S{x — tk'), are not allowed since they are unbounded, even though they obviously satisfy 
{d{tk/, z), f) — f{tk') (they also cannot be considered reproducing kernels for the same reason). 

Notice that a GLLSQ setting and its solution epitomized by (IA-4P entails little more than a positive 
definite norm and inner product structure. Both DIDACKS and RKHS approaches simultaneously display 
both GLLSQ and standard collocation aspects; hence, there is an underlying Gram matrix in both cases. 
This dual aspect of both DIDACKS and RKHS approaches comes together in Definition A.l, but there a 
number of side issues to be dealt with. For example, while a choice of orthonormal basis functions usually 
makes sense for many GLLSQ processes and the trivial situation T = I occurs in this case, for RKHS 
and DIDACKS approaches matters are somewhat different. Although a set of orthonormal basis functions 
may be known for the family of admissible functions and it may be used to reexpress the replicating or 
reproducing kernel, as discussed below, need not necessarily exist in such cases and the actual cases 
where T = I occurs are very much the exception rather than the rule. One case where T = I could occur in a 
reproducing-kernel like environment is if K{z, t) were a Dirac delta function, but as just noted a Dirac delta 
functions are not an allowed kernel choice. Likewise, in the DIDACKS environment a set of point sources 
located on dH. would produce a diagonal T matrix, but in DIDACKS theory it is assumed that the kernel 
functions are harmonic in not only the interior of Q, but also on dfl so this is not an allowed case either. 

Before proceeding, notice that there is an entirely different standard strategy for forcing a candidate 
linear combination of basis functions (p to match some given function /, which is characterized by directly 
imposing the function matching condition implied by (|A-ip . psp or some analogous relationships and then 
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simply solving the resulting linear equation set for the coefficients of Lp. In general, no inner product structure 
is required to perform this standard type of collocation fit, so there is usually no associated Gram matrix. 
In the literature such approaches are referred to by the broad rubric collocation, but here is it useful to 
label such approaches direct collocation, when the need to distinguish them arises at all. Other approaches 
also use this same strategy, but do not use the label collocation. One such approach is MFS, as commonly 
understood, which was addressed in Section [6] and where connections to material in other sections of this 
appendix were pointed out. 

To put all of this in proper context, first an analysis and extension of (jA-ip is called for, in order to better 
correlate it with Definition A.l. Since a general context is still being assumed, it is useful to first break out 
the collocation function matching condition associated with (jA-l[) and introduce a separate name for it. 
Thus the criteria that some arbritrary, but admissible, function g(i) with i G $7 matches some prescribed 
constants at the iV^ specified points 

f{tk) = Ck , (A-8) 

is labeled the point collocation condition. For purposes of discussion only one condition of the form (|A-8P will 
be entertained at a time so that the Ck are to be regarded as fixed real constants throughout the discussion. 

For 77 e 57, let the values of 77 that satisfy f{ri) = be called node points. Obviously, when two separate 
functions / and g both satisfy the point collocation condition g{tk) — f{tk) and one of these functions will 
be said to satisfy the point collocation condition with respect to the other. In this case / and g can differ 
by an arbritrary admissible function x(Oi 5 = / + Xi provided that the tk are node points of x- x(^fe) = 0- 

Next, in order to amplify on the foregoing and to consider it from a slightly different perspective, first 
notice that when U = U{z) is expanded as a linear combination of Nk basis functions the coefficients are 
uniquely determined by the point collocation condition provided that the resulting system of equations is 
invertible. Here this collocation fitting process is often considered without assuming explicit norm and inner 
product structures. Specifically, let {/ifej^i be a set of suitable basis functions that are not necessarily 
orthogonal. The term suitable here generally means that this set of basis functions are linearly independent 
and are appropriate to the situation being modeled. Then, as in (TO]) , let 

U{z) = Y,^lkhk{z) . (A-9) 
fe=l 

(|A-8|) then becomes 

Ha< = C (A-10) 

where H is the matrix whose components are Hkk' — hk{tk') and C is the vector whose components are 
Ck- For the sake of discussion assume always exists. Here the number of basis functions used has 

been constrained to be the same as the number of sample points on the RHS. In general, if the number of 
sample points is larger than the number of basis points then a LLSQ fitting process is implied rather than 
the process (jA-10|) . while if the number of basis functions is larger than the number of sample points then 
the associated H matrix does not have full rank. Both of these cases are of no particular interest here and 
are not considered in the sequel. If the hk are expressed as is, and not in kernel form, (|A-9[) will be called 
a direct collocation process. It also represents what might also be labeled a forced least squares LLSQ fit, 
where the number of sampled residuals is Nk- Next suppose that a kernel, 7i(/9, z) can be associated with 
the /ifc's by the ansatz hk{z) = Ti-ipk, z). The idea here is that when a proper kernel fmiction exists the 
index k can be naturally transferred to a set of indexed points {pk}k=i^ where it is assumed that pk ^ Pk' 
whenever k 7^ k' . Here pk G ^ need not be assumed, but pk € K" generally is assumed. In this case and 
ignoring ill-conditioning issues, since Hkk' — 'H{pk, tk') the main criteria for a successful kernalization of 
hk is that H"'^ exist. The matrices H and T obviously play the same role, so only the matrix T will be 
referred to in the sequel. As noted above, when neither an associated reproducing or replicating property nor 
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an associated norm is assumed, then (|A-9I) exemplifies a direct collocation kernel based process and direct 
collocation, in one form or another, is pervasive in the applied scientific literature where it is usually simply 
called collocation. 

Reproducing kernel and replicating kernels fits can be considered kernelized examples of (|A-10[) . First, 
as an aside, evaluate Equation (|A-6|) at the points tk> to obtain 

fe=l 

This equation, in conjunction with Equation (|A-7p . then implies that (p satisfies the point collocation condi- 
tion with respect to /, since K{tk, tfc') — K{tk', tk) holds by assumption. 

There is one additional distinction to be considered for replicating and reproducing kernels with regards 
to the constraint condition (jA-Sp . As noted, for these cases, while it is a Gram matrix, the T matrix is only 
indirectly related to the set of orthogonal basis functions that span the family of functions and, moreover, is 
always a square matrix of size Nk for any specified set of Nk evaluation points. If only the general form of the 
kernel function is known, then unless the explicit dimension, Nf, of the span of the family of basis functions 
is known then > Nf may well occur. In this case (jA-81) either completely determines or overdetermines 
the coefficients of ip, which is a misapplication of the underlying technique. In this situation (jA-Sp will be 
said to overwhelm the family of admissible functions. A separate condition must be imposed in this case since 
it is a property of the class of admissible functions and not necessarily the kernel basis functions themselves, 
but it is dependent on the number, Nj^, of evaluation points. Crudely put, as discussed in Section U] and 
below, the requirement is that Nf > Nk- (It will soon be painfully apparent that the real relevant concepts 
are that of uniformly N pointwise independent family of admissible functions over Q, which was introduced 
by Definition m 1 , and it is difficult to avoid it.) 

Much of the foregoing discussion can be understood by examining a reproducing kernel collocation fit 
and in the ensuing discussion it should be plain why certain conditions (such as |T| 7^ 0) must be assumed. 
Suppose, as usual that is the region of interest. First, consider the kernel form: 

i^(z, i)=X^^j(^)u,(i) , (A-12) 

where {uk}f^i are a complete set of orthonormal basis functions over fl: 

{Uj, Uy)^ = Sjj' , (A-13) 

where as in [T], the bound variable(s) for the inner product in question are indicated by a subscript. Notice 
that K obviously has the reproducing property: 

{K{z, t), /(z)), = Y.^,{t) {u,{z), f{z)U = fit) . (A-14) 
fc=i 

A reproducing kernel expression of this type will be called regular. 

An RKHS fit corresponding to (jA-6p and (jA-7p can be rewritten in a revealing way for regular kernel 
expressions. First rewrite (|A-6I) as 

k=l k=l 
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where 

Vk{z) = K{z, tk) = y^^Uj{tk)uj{z) 
i=i 

or in matrix form as 

V = WU (A-16) 

where V(z) is the vector whose fc'th component is V}^{z)^ U(z) is the Nu vector whose fih component 
is Uj{z) and W is the N]^ by Nu matrix whose components are Uj{tk)- Equation (|A-16[) clearly represents 
a transformation to a new set of basis functions. Unless W, which does not depend on z, has rank the 
new set of basis functions will not be linearly independent. Observe that (|A-7|) can be written in the form 

J^Gk'k^ik^ f{tk') (A-17) 

k=l 

where Gk' k {vw, Vk)z corresponds to a Gram matrix G and this matrix will be non-singular only if W 
has rank Nk- Further notice that the rank of W may depend on the actual interpolation point locations tk, 
which corresponds to the "geographic factors" discussed below. 

Point Collocation Conditions With Operators 

Next consider the analog of (jA-Sp for operators that is implied by the point replication condition Definition 
A.l. Thus suppose that this condition can represented by the action of a global linear differential operator 
Lk that may differ from one k value to another and that it can be expressed as a sum of partials of various 
orders with respect to the components of t at each location. Within the context of Section [51 this means that 
if £;fc is the operator associated with k then it can be represented as a sum of partials of various orders with 
respect to the components of Pk- For concreteness, consider the case where a single differential operator 
exists for all k values so that iLt = Hk and 

^fe/ ^tf{t) = Ltfitk) ■ (A-18) 
t=tk 

By extension, the notation for the situation where Lt varies from one set of fc's to the next will be denoted, 
for convenience, simply as 

^kf{t) = Lkf . (A-19) 
t=tfc 

The point replication analog of (jA-8|l then takes the form 

= afc, , (A-20) 

for fixed real constants at- Just as before, when g satisfies the condition l|A-20[) . it can be written as g = Z + x 
where the tkS are now node points of ^k{x) rather than x- Expression (|A-20p will be called the generalized 
point collocation condition. 

At the risk of being overly abstract, the exact nature of the linear differential operator and its 
relationship to the condition of Definition A.l has been temporarily left open to allow for the greatest 
generality. Presently any residual confusion that enters on this account will be dealt with. For now can 
be taken to be an undefined primitive or it can be considered to be a multiplicative constant, whichever the 
readers wishes. 
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Since the conditions implied by Definition A.l hold, it follows that 



{Bk, f) ^ iikf = . (A-21) 

Also (|A-3p and (|A-5p hold. Here (p obviously satisfies the condition implied by (|A-20p since (|A-5I) . as 
expressed in the form (j A-4[) . can be immediately be rewritten as {Bk', — {Bk', /) so that condition (|A-20[) 
expressed for /, along with (|A-2ip . also implies it holds for ip: 

i^fct^ = ak . (A-22) 

Likewise multiplying (jA-Sp . as expressed in the form (jA-4p . by ^k' and summing over k' gives 

^) = {p, /) , (A-23) 

which can be used to eliminate {tp, /) on the RHS of (|A-2p : 

\\f^V'f = \\ff-M'- (A-24) 
Since the norm has been assumed to positive definite obviously 

11/11 >ll¥'ll- (A-25) 

It has been assumed that the ak and / and thus ip are fixed. Consider some other arbritrary admissible 
function that also satisfies condition (|A-20p : 

Hkg = ak , (A-26) 

then if ||g|| < ||/|| is it possible for g to be less than ip7 No, for the following reason. From assumptions 
about the nature of the action of the operator and inner product replication condition (Bk, g) = ^kg, so 
that (|A-26p implies g satisfies all of the conditions assumed for / starting with (jA-20p and thus that / can 
be replaced by g up through and including (|A-25P : 

Il5ll>ll^l|. (A-27) 

There are several side issues that need to be considered; however, most of them will be addressed after the 
main theorems are obtained. 

An attempt will now be made to summarize the progress so far in several definitions and theorems. While 
the requirement that the underlying vector space of admissible functions have at least dimension Nk might 
seem to be sufficient, as previously indicated, subsequent discussions and examples show that it is not and 
the stronger condition that the family of functions is uniformly Nk pointwise independent will have to be 
assumed. 

Before proceeding it is useful to clarity how the linear differential operator Lk enters here. Towards that 
end some additional precision in terminology is in order: 

Definition A. 2 A linear differential operator Lk{t) that is associated with Definition A.l as discussed 
above, will be called k-universal when it is the same basic operator for all fc, which is to say it has the same 
form for all values of k [where the evaluation convention implied by (jA-19p obviously is still assumed to 
hold]. 

Next, a definition is introduced to limit the scope of operators under current consideration: 

Definition A. 3 A linear differential operator ilfc(i) will be called primitive when it is fc-universal and it has 
either the form Lt — X ov Lt = A|i|, for some non-zero real constant A and < |t| < oo. From the evaluation 
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convention mentioned in Definition A. 2, this means that for A 7^ either — X for all k or Lk — X\tk\ for 
all k and < \tk\ < 00. 

In Definition A.l, while an explicit side mention is made of the point replication property satisfied by 
all DIDACKs, it clearly also applies to reproducing kernels. When Definition A.l is understood in the 
sense of a primitive operator, as outlined in Definitions A. 2 and A. 3, the various cases of immediate interest 
for standard RKHS and scalar point source fits can be dealt with. For example, a RKHS collocation fit 
corresponds to ZLj. — 1 for all k and p7p implies that a scalar point source fit corresponds to L^. — |tj.| in 
the current notation. It is worth noting that if a linear differential operator is not suitable, then there is 
a real chance of obtaining meaningless and trivial results. For example, if Lkf 0, then (jA-27p becomes 
\\g\\ > 0. The index k dependence of Hk also leaves open "geographic implications" that must be considered. 
("Geographic implications," are discussed below in a related context.) These geographic implications can 
be overcome by generalizing the idea of uniformly N pointwise independent families of functions to include 
the action of operators. For example, if it can be shown that for some /c-universal operator H, Lf maps a 
given family of uniformly pointwise independent functions into a corresponding family of Nk pointwise 
independent functions, then the resulting T matrix will be invertible and the replication or reproducing 
property (jA-22p will be non-trivial. 

For general linear differential operators that are non fc-universal issues that must also be considered. 
Observe, for example, that L^. itself may be a linear combination of other differential operators and thus 
it may have internal parameters. If these other parameters are to be fit, then it is a simple matter to 
separate out the resulting relevant types of operator basis functions, but in this case having multiple types of 
operators associated with a single location produces non-trivial complications because it is generally harder 
to prove that these various types of operators all act independently of each other. An example of one such 
situation arises is for a combined point mass/dipole associated with location Pk = Xk/r\. It is obviously 
much harder to prove the linear independence all of the component parts of this sort of combination than it 
is when only simple operators are present; moreover, it unclear how to even generalize existing point mass 
linear independence proofs when even only simple dipole terms appear. (Current plans are to present a 
more general line of proof that includes dipole linear independence as a special case in a future article.) 
Other complexities also enter when considering compound linear differential operators, since they may well 
not actually yield linearly independent basis functions and, in fact, the issue of linear independence may 
depend on the underlying properties of the family of admissible functions itself. For example, first let 
(Di, ©2, 'D3)"'" := {d /dx, d /dy, d /dz)^ — V, then consider the operators Lij := DiDj. There are at 
most six independent components since (for any well behaved family of functions) Lj i — Lij\ however, if 
/ is harmonic then there will be at most five independent operators of this type at any given point due to 
Laplace's equation. For this example, the interaction of the operator form Lij with the class of admissible 
functions is fairly obvious, but for certain operators and classes of admissible functions this interaction may 
not be so transparent. (For example, as pointed out later Dr/ is not harmonic when / is harmonic and 
thus it is not in the family of admissible harmonic functions defined over, say, Q,q C R'^; hence, it cannot 
rigorously be represented as a harmonic expansion.) 

For the above reasons, as well as ones discussed below, it is expedient to restrict ourselves to relatively 
easily obtained results. When this is done, all of the above can then be summarized in four main theorems. 
Recapping the line of argument leading up to (|A-22p then yields: 

Theorem 3. When the generalized point collocation condition is satisfied by a primitive linear operator and 
the conditions of Definition A . 1 are also satisfied for a uniformly Nk pointwise independent kernel environ- 
ment, then the resulting fit, Lp, also satisfies the same generalized point collocation condition; moreover, a 
point collocation condition of the form ^ A- 20(1 is also implied. 
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This theorem guarantees that the given data set is repKcated by the resulting fit. Here, notice that if 7^ 
then for primitive operators ak = ^f{tk) or ak = X\tk\f{tk) in (jA-20|) and thus the condition (jA-2ip imphes 
the point coUocation matching condition (|A-8p and conversely. 
The line of argument leading up to (jA-27|l yields: 

Theorem 4. As in Theorem\^ given that the generalized point collocation condition is satisfied by a prim- 
itive linear operator and the conditions of Definition A.l are also satisfied for a uniformly Nk pointwise 
independent kernel environment, then the resulting fit ip also satisfies the \\g\\ > \\lp\\ for any other admissible 
function g that satisfies either the same generalized point collocation condition or point collocation condition. 
(Equivalence of the conditions iA-S^} and llA-20\) is also implied in this case). 

Next it is useful to consider Theorems [3] and [4] specifically in the context of DIDACKS point mass fits: 

Theorem 5. R'^ DIDACKS point mass fits replicate the specified data points of f used; moreover, the 
resulting fits minimize not only \\f — y^\\ but \\ip\\ is less than or equal to that of any other admissible function 
that also matches the .specified data points. 

Proof. As previously noted in Section |4l there is a proof in the literature that a finite set of point mass 
basis functions are linearly independent for the situation where the sources are contained inside some sphere 
and the harmonic region consists of the exterior of this sphere. If linear independence of point mass basis 
function can be shown in the complimentary case where the sources are in VIq and f2[o] is the region of 
interest then all of the conditions for Theorems [3] and |4] to hold will obviously be satisfied and thus the 
desired result will follow. Towards this end consider the method of images for the R3 sphere [U. This 
technique states that to each charge in the interior of the sphere there corresponds a charge of a prescribed 
strength of opposite sign and at a determined location such that the net effect of these two changes together 
gives a zero potential on the surface of the sphere in question. [The relative positions of these two charges is 
prescribed by PH)) .] Suppose to the contrary that a counter example did exist where a collection of non-zero 
changes in r2[o] produced a zero potential over 9r2[o]. Then from the method of images each individual charge 
in this combination would have a image charge that also produces a zero field over 9f2[o] and, furthermore, 
the sum of all of these charges would produce a linearly dependent counter example to the existing theorem 
for point mass basis functions. From this contradiction the desired linear independence result follows. □ 

Geometric Interpretation of Collocation Fits 



Here, temporarily leaving aside the issue of general operators, it is worth observing that it is possible to 
give an explicit geometric interpretation for collocation fits, but the correct one appears to be slightly more 
involved than the generally well-known one. Before proceeding, as a side issue, first consider the geometry 
that occurs in n-dimensional Euclidian space (K") when the values of iV < n + 1 vectors are specified as 
a constraint condition and one wants to find the vector of minimum length that satisfies this constraint. 
(As in other parts of the article, at several places here the reader may wish to draw a diagram for her or 
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himself.) Since A'" points determine an iV — 1 dimensional hyperplane, this condition of constraint determines 
an A'^ — 1 dimensional hyperplane. The vector of minimum length from the origin to this hyperplane is then 
the one that is orthogonal to it. There is, however, a notable difference between this situation and the 
present one exemplified by (|A-8|) or (jA-lOp — equations (|A-8p and (jA-lOp are conditions of constraint that 
directly determine the various basis function coefficients themselves. Thus the closest ordinary vector analog 
of (|A-8p is the one where the components of some arbritrary vector are fixed with respect to a subspace 
spanned by a set of basis vectors (corresponding to the various chosen basis vectors, hk, themselves), but 
not fixed with respect to the whole space. To see this, for example, consider the following hypothetical 
situation that uses basis functions that have a coordinate representation. Thus suppose that, contrary to 
the assumption of boundedness, Dirac delta functions are allowed basis functions and that an n-dimensional 
volume integral of f g specifies the inner product for / and g, then the minimum norm solution to (jA-8p 
that results is obviously f{z) = C^kS{z — In this degenerate case the 5{z — ffe)'s are the basis 

functions and the C^'s are the projections of / onto them. Thus (|A-8p is not a condition that specifies a 
hyper-plane of dimension A^^ — 1 in the function space, but rather it is a condition that implies the projection 
of / onto a set of basis functions, which may be as yet undetermined. For this whole geometric situation 
not to collapse, it is necessary that Nk < n hold. Here Nk = n is not of interest since it completely fixes / 
and g = f — (fi. To see that the proper dimension is indeed only Nk and not Nk — f , in (|A-9p and (|A-f Op 
consider a function of the form U (z) = X^^i l^k5{z — tk). In this case H — I and each individual constraint 
condition in (|A-10p fixes the magnitude of the associated /ifc for the basis function 5{z — tk). Hence (lA-IOp 
is not a condition specifying the points themselves (i.e., vectors) in the function or Hilbert space, but it is a 
condition specifying the components of the basis vectors. 

Finally, observe that if additional Dirac delta basis functions of the form 5{z — Zk') for Nk' > k' > Nk are 
added, where Zk' 7^ tk for all k' and k, then the associated coefficients fik' are not in any way constrained. 
Specifically, let U'{z) denote this new functional form with the added basis functions, then at any of the 
evaluations points of (|A-10p . say p: 



Only the components of U{z) that fall within the constraint subspace are fixed and the other projections 
of U (z) into the compliment of this space is not fixed at all — some sort of least squares or minimum norm 
solution is needed to fix that part of U (z) that is in this complimentary space. Clearly, in so far as possible, 
these minimum norm solutions will tend to minimize the part of U{z) that is in this complimentary subspace 
so that the chosen solution will be the one where ^k' = for all Nk' > k' > Nk- Within the context of 
equation (jA-27p . U{z) (p and U'{z) g{z). Observe that when overlapping basis functions — such as 
fundamental solutions — are used this subspace decomposition is obviously not nearly as clean. 

With the foregoing as backdrop, consider now the situation at hand. First, take two different arbritrary 
admissible functions / and g that satisfy (|A-8p . Here x, the difference of these two functions, was introduced 
above. When either f or g is fixed and the other is varied, X = f ~ 9 spans a linear vector space since 
aXi(ifc) + bx2{tk) = for a and 6 G M, provided of corse that Xi{tk) = Xi(^fc) — 0- Let {x} denote the span 
of this full vector subspace. Temporarily leaving aside the question of uniqueness, of all possible functions 
that satisfy this condition one will have minimum norm; let ip denote this function. The requirement that ip 
yields the minimum \\(p\\ of all those functions that satisfy condition (|A-8p is equivalent to the requirement 
that it must be orthogonal to the any vector in {x}'- (Vi x) = or (ip, /) = {p, g). Suppose / is the fixed 
element in the representation of {x} and g is the variable one. Clearly g = ip is one possible choice so 




k=l 



k'=Nk + l 




(A-28) 
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(^5, Lp) = (ip, /) must hold, which is (|A-23|) . It is then obvious that ||(^|| is a minimum from the steps leading 
up to jX^ . 

When the CFGLLSQ condition holds for primitive operators, solutions that satisfy the collocation con- 
dition also are GLLSQ fits. Suppose that / is to be fit with a linear combination of the form (jA-Sp . Thus 
11/ — ifW^ is a minimized when the projection of / onto the subspace spanned by the set of basis functions 
{'^fc}£i is matched exactly in the sense of (|A-8|) . Geometrically this occurs when f — ip is orthogonal to 
the plane spanned by the set of basis functions {Bk}^!li- {f, = (<p, /)• Furthermore, since f — p is 
orthogonal to this subspace it must also be orthogonal to each of the component "basis vectors" of this 
subspace: {Bk, f — f), which imphes (|A-4p . 

A Simple Example and Related Issues 



A brief discussion of counter-examples and associated issues is in order so as to make it plain why certain 
conditions (such as |T| ^ 0) must be assumed (if it is not already obvious enough to the reader). Setting 
the stage just as before, first consider the kernel form: 

K{z,t)^J2uk{z)uk{t) , (A-29) 
fc=i 

where < Nf + 1 and {uk}^!^ are a set of orthonormal basis functions: 

(wfc, Uk') = 5kk' ■ (A-30) 
When Nu = Nf , K is called regular and obviously it has the reproducing property: 

Nf 

(Kiz, t), f{z)), = W (^fc(^)' /(^))- = 

fc=i 

where as in [1], the bound variable(s) for the inner product in question arc indicated by a subscript. If 
Nu < Nf a projection onto the subspace of {/} may occur where {/} denotes the complete linear span of 
the family of admissible functions. In this case K will be called a pseudo-reproducing kernel and such a 
kernel may be mistaken for a full reproducing kernel unless careful attention is paid to make sure that the 
kernel basis functions and / are fully representative members of the family of admissible functions. This 
may happen, for example, when K is composed of basis functions obtained through the action of a linear 
differential operator. Such an operator may also map {/} into some other span of functions: {£jf} ^ {/}. 
An associated undesirable situation occurs when / is restricted to be in some sub-span of the family of 
admissible functions, say {wfej^jy^ where A'^i > 1 or < Nf. Again, while the inner product of a function 
of this type and a regular reproducing kernel yields a function of this type, it may also be the case that 
a pseudo-reproducing kernel gives the same result and that simply showing that the reproducing property 
holds in this situation may not be sufficient. For these reasons it is clear that the reproducing kernels and 
/ must be general members of the family of admissible functions. 

One concrete elementary example should suffice here. For simplicity only one-dimensional examples need 
be considered, so let x, y € M represent the two independent associated variables over Q, rather than z and t. 
Further, let T„p denote the family of polynomials in one variable of degree less than or equal to Up for some 
fixed interval, say [—1, 1]. Since the set of Legendre polynomials, P„(a;), are a well-known set of orthogonal 
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polynomials over this same interval with inner product specified by 

1 

(/, g):= J f{x)g{x)dx, (A-32) 

-1 

they can be used to construct the required orthonormal basis vectors and associated regular reproducing 
kernels: 

/2fc — 1 

Ufe(x) Pfe_i(a;)y . (A-33) 

First consider a pseudo-replicating kernel constructed with this set of basis functions via (|A-29[) with Nu = 2 
and let Nk — 3. Since Pq = 1 and Pi = x, ui = l/\/2 and U2 = a;-\/2/3, Tkk' = 1/2 + {2/3)xkXk' ■ Consider 
the sample points xi — 0, X2 = — \/3/2 and 2:3 = \/3/2. It is easy to show that |T| = for this case. This 
example leads to other questions about pseudo-reproducing kernels and over-represented cases in general, 
so these cases are best simply avoided at the outset. Clearly the linear differential operator d /dx reduces 
the dimension of the linear span of by one and, when applied repeatedly, will lead to over-represented 
cases (i.e., cases where > Nj) for any fixed set of polynomials. Observe that this same concern occurs 
with regards to a harmonic family spanned by spherical harmonics of finite degree and order over 17 [0] since 
such spherical harmonics can be reexpressed as polynomial expressions involving finite powers of x, y and 
z. For the infinite dimensional case the action of a linear differential operator may also reduce or alter the 
span of the family of admissible functions. Any proof of the collocation and minimum norm property that 
omits the stipulation that the underlying linear system is invertible and not over represented is immediately 
suspect — especially, if it purports to show the minimum norm property for linear operators. 

It is clear that questions about the nature of the T matrix arise when several of the tk are node points 
or especially when all, or many, of the basis functions share some common node point, since all of the 
associated basis function terms drop out of the T matrix. (The converse here is also interesting since |T| 7^ 
implies limits on possible shared node points of orthogonal basis function sets.) Given that only the form of 
the kernel function may be known and that connections to the form (IA-29P may be unclear, especially for 
DACKs, this whole set of issues is best avoided. For all of these reasons, simply requiring that exists 
and that Nk < Nf is much easier than explaining away all the special cases, even when it can be done. 



Geographic and Other Factors 



There are also other factors here to consider that have not, as yet, been properly addressed. Since the 
primary setting of this article involves M.^ harmonic functions, some of these additional points will be made 
specifically with regards to harmonic functions. First it might be thought that simply requiring harmonicity 
is, of itself, sufficient, but at the very least, it should be stipulated that the region of interest is contained in 
M" for n > 1, since harmonic solutions in M} satisfy d^ f /dx"^ = 0, which has a general solution that spans 
only a two dimensional space and so n/ = 2 (i.e., f{x) — ax + b for the case x G [c, d]). For the rest of this 
discussion only the M'^ harmonic case will be the primary focus. Second, overall constant factors are often of 
particular concern, especially if linear differential operators are involved. For example, suppose that as for the 
I norm, harmonic functions over fig that fall off sufficiently fast as r ^ 00 comprises the class of admissible 
functions, then in general Zf may not be in this class of admissible functions itself — especially if it is an 
operator that involves integration. Third, V'^ f — does not imply W^d f/dr — since d f/dr = r~^X -Wf 
(which follows form basic vector analysis or an application of the chain rule) Fourth, saying that / = 

*Here an interesting example arises from the spherical approximation of gravity anomaly Ag, which frequently occurs in 
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over Q implies certain geographic or coordinate dependent considerations. These "geographic issues" are 
worth considering separately in some detail. 

Two entirely different "geographic factors" will be touched on here. The first concern is a general one that 
is somewhat independent of harmonic function considerations. The concern here is that reproducing kernel 
counter-examples can be cobbled together for admissible families of functions that are separately defined 
over individual partitions or subregions of fl, such that each subregion has a different Uf; consequently, 
simply saying that Nk < Uf is not sufficient unless Uf is taken to be the smaller of these two values. As an 
example, let the family of admissible functions be defined as polynomials of fixed degree and order in x, y 
and z for the interior of a unit cube or sphere, such that /(x, y, —z) = —f{x, y, z). Suppose that the highest 
allowed order is different for z > than it is for z < 0. If a reproducing kernel can be defined for 2 > 
that is different than it is for z < and if Uf depends on whether the sample points Pk are located in the 
upper or lower half-plane then a counter-example can be constructed that displays geographic dependence. 
It is unclear whether some sort of family of admissible functions and kernels exists that displays a similar 
geographic variation that is continuous and gradual. 

The second geographic factors center on issues related to (N) pointwise independence, but it is useful to 
first consider the issues involved independently form this concept in oder to see that in the end for expediency 
that this concept (or something like it) is probably required. From the uniqueness of Dirichlet boundary 
conditions, if / is specified over some closed surface that is contained inside f2, then it is fixed in the interior 
region bounded by this surface; moreover, to the extent that harmonic extrapolation is well defined, its 
values are also fixed throughout the whole of Q. (This issue of harmonic extrapolation is taken up in greater 
detail in Appendix B.) This clearly also implies that ip is, to some degree or another, fixed by one form of 
DIDACKS replication result or another and it brings up the possibility of whether there is an associated 
geographic constraint. The question as to whether harmonic functions are geographically constrained or not 
may, at first, seem superficial, but it is not immediately apparent that there is a sufficient amount of implied 
freedom in the concept of an M*^ harmonic function to guarantee that a harmonic function / in 17 always 
exists such that each value of f{Pk) can be assigned an independently for the set of arbitrarily specified Nk 
distinct points Pk G It is assumed that < \Pk\ < 00 regardless of fl. If this property holds for harmonic 
functions defined over ri[o]j then by the Kelvin inversion theorem it must hold for harmonic functions defined 
over rio as well. The same can be said for fli, since the source regions are assumed to be bounded in this 
case as well. 

This property does, in fact, have some surprising implications so it is not trivial. In particular, consider 
what this implies for harmonic functions defined over Qq. First, in analogy with Tabled let a S/e gridding 
for a sphere of radius Rj be a set of points on this sphere that are packed as closely as possible with the 
requirement that the closest two distinct points have an angular separation of 6. Then e is the resulting 
maximum angular separation of any two adjacent points and by closest packing it is meant that e is as small 
as possible, in some sense or another. Practical examples of i5/e-like griddings are given in Table [21 Clearly 
there is a finite number of points in a 6/e griddings for each specified S. If admissible families of harmonic 
functions are uniformly pointwise independent then the following theorem obviously holds: 

Theorem 6. Let S denote an arbritrary, but fixed, positive parameter and let {Rj}^J^i be a set of Nj real 
numbers such that 00 > Rj+i > Rj > 1. In M.^ , given a S/e gridding over each sphere of radius Rj, a 
harmonic function, f, in Qq exists that will match each of these S/e gridding points when arbritrary values 
are assigned to f at each of these points. 

physical geodesy. In this approximation Ag := dV/dr — {2/r)V where V^V = 0. Although, ^^dV/dr ^ 0, using the product 
rule for the Laplacian it is easy to show that W^dV/dr = V^[(2/r)V] so that V^Ag = 0; however, some care is required here 
since, for example, V^[dV/dr - i2/R)V] = S7^[dV/dr ^ 0. 
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This theorem is somewhat surprising, since one might reasonably expect from the uniqueness of the 
Dirichlet boundary value problem that specifying values on a dense grid at a sphere of radius i?i > 1 is 
sufficient to fix the values on all other surfaces; however, Theorem [5] implies that when some fixed 0, (j)) 

is given and it is sampled at the points of a (5/e grid at a radius i?i > f specified by an associated 5, which 
can be any small fixed number, then other harmonic functions necessarily exist that match / at the given 
points Ri and take on arbritrary values at specified points on a larger sphere (or even on an array of 
spheres). Moreover, from the results discussed above, ^ provides a constructive way of obtaining a function 
of minimum weighted energy norm that will match all of these prescribed points. As a specific example, 
consider a fit over various concentric spheres to the nonharmonic function / — r™ with m > 0. Here as 
(5 ^ it is easy to show that although the collocation (or interpolation) values match at collocation points 
to machine accuracy, away from the collocation points yields a worse and worse approximation as (5 — > 0, 
as one would expect, since the natural attenuation of a harmonic function is obviously not being modeled. 
(As elsewhere, it is instructive and fairly easy to implement this sort of example and readers are encouraged 
to test this sort of thing for themselves.) 

From all of the foregoing analysis it is apparent that the best way to make progress is to first settle on 
some specific allowed attributes that the operators and family of admissible functions must posses and then 
state theorems that can be easily proven for operators and admissible families of functions with these specific 
traits. From all of the issues brought up above, it would seem to be apparent that it necessary to require 
that the family of admissible functions be either uniformly pointwise independent or uniformly N pointwise 
independent, in accord with Definitional. Future progress in this area beyond the results presented here 
can then be undertaken by introducing new definitions for linear operators, as the need arises, that are 
generalizations or distinctions, in some sense or other, of those operators characterized by Definition A. 4, 
given below. The above discussion and counter-examples should then serve as a backdrop for these possible 
future considerations. 

Formal Result 

As just noted, in order to introduce a certain measure of clarity and closure it is expedient to specify 
concrete traits of not only the family of admissible functions to be considered, but of the allowable operators 
for point reproducing and replication conditions. Since relevant results can easily be adapted to the case 
of uniformly N pointwise independent families of admissible functions, if they hold for uniformly pointwise 
independent functions, only uniformly pointwise independent functions will be considered in the rest of this 
appendix. 

With regards to the inner product action of point reproducing or replication basis functions, a slight gen- 
eralization of primitive operators will be considered, which is labeled semi-primitive. Thus for an appropriate 
inner product environment, consider basis functions Bk{X) that have a point reproducing or replication prop- 
erty with an associated semi-primitive operator, so that the point reproducing or replicating property can 
be written, in analogy with (|A-21[) . as 

{Bk, f) = h{Pk)f{Pk) , (A-34) 

where h{Pk) ^ 0. (Here for a primitive operator, h(P) has the form h{P) = A or /i(P) — A|P| for some 
nonzero A.) Formally: 

Definition A. 4 A linear operon L that is associated with some family of admissible functions "3 defined 
over some subregion 57 of R'^ is said to be semi-primitive if for all X G 51, L, — h{X) or Lf{X) = h{X)f{X), 
where h{X) ^ 0. Likewise, a set of basis functions {Bk}^^i in 3^ are said to be semi-primitive if they can be 
characterized by the action of a semi-primitive operator, as in (|A-34[) . 
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Here, for concreteness, 3^ was defined over some subregion of M'^, but the same definition also holds for 3" 
defined over any suitable R" domain. 

Since 3^ is assumed to be uniformly pointwise independent, the Bk specified by Definition A. 4 must be 
linearly independent. Suppose to the contrary that the Bk's are linearly dependent and thus that some set of 
gfc's exist such that J2k=i ^k{X) = 0, where at least one of the q^'s are non-zero. Specifically, let q^' ^ 
and then reorder the basis functions so that Bk' is the first one and thus qi ^ 0. Expression (|A-34p implies 
then that Y.kli<lk{^k, f) - Ekli 1kh{Pk) f (Pk) = 0. Here choose / such that f{Pk) = i if fc - 1 and 
f{Pk) = 0, if /c > 1, then X^^i <lkh{Pk) f (Pk) = qih{Pi) ^ 0, which is a contradiction. Using a more general 
notation, in summary: 

Theorem 7. // the family of admissible functions 7 is uniformly pointwise independent and the set of basis 
functions {Bk}'^^i are semi- primitive reproducing or replicating kernels, which is to say they have an inner- 
product action characterized by !lA-34^ , then this set of basis functions {Bk}k^i is linearly independent and 
the resulting T matrix is invertible. Conversely, if for any finite Nk a semi-primitive basis function inner- 
product expression of the form llA-34\ l holds for {Bk}k=i o.^T'd if the resulting T matrix is invertible for any 
set of distinct points {Pk}k=i located in the domain of interest, then 'J is uniformly pointwise independent. 

The proof of the first part of this theorem was just outlined above and its converse, which is the second part, 
is obvious. 

Thus if it can be shown that harmonic functions are uniformly pointwise independent over some DIDACKS 
region of interest, say ri[o], then it follows that an associated set of scalar point source basis functions must 
be linearly independent. As previously indicated, a proof exists in the literature that scalar point source 
basis functions are linearly independent; however, another proof of this fact would be enlightening so an 
attempt is made in Appendix B to show directly that harmonic functions are pointwise independent. While 
not entirely successful, this attempt reveals certain hidden dangers in the process of harmonic extrapolation 
(i.e., downward continuation for unbounded exterior cases or outward continuation for bounded interior 
cases). 
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Appendix B 



Attempt to Prove Harmonic Functions are Pointwise Independent 

This appendix makes an attempt to directly prove that harmonic functions defined over an appropriate 
region Q are pointwise independent. The DIDACKS point repUcation conditions for a given region Q, in 
conjunction with the fact that point source basis functions are Unearly independent (which, as previously 
indicated, has been proven in the literature) imply, by Theorem [7] at the end of Appendix A, that harmonic 
functions are pointwise independent. Although the pointwise independence of harmonic functions is not in 
question, as noted elsewhere in the article, it is a very significant attribute of harmonic functions and thus 
one might wish for, and expect to find, a direct proof of it; consequently, this appendix gives an attempt to 
supply such a proof. Also, observe that the DIDACKS point replication conditions for a given region fl, in 
conjuction with Theorem[7]at the end of Appendix A, imply that if the pointwise independence of harmonic 
function over Q can be shown, then this constitutes an independent proof that scalar point source basis 
functions are linearly independent, which would also clearly be of interest. Various side issues and auxiliary 
theorems of note will emerge as the attempt to prove that R'^ harmonic functions are uniformly pointwise 
independent unfolds. 

This attempt will be approached by way of several lemmas and side theorems specifically for the il[Q] 
domain. As previously noted, the results obtained for this special case can be extended to other domains 
of interest. For convenience, let 7i[o] denote the family of associated harmonic functions. Recall that these 
functions are bounded and harmonic on the boundary dQ[o] well as in the interior; moreover, it is assumed 
that they have bounded partials [they are also smooth (i.e., C°°)]. Since a harmonic function is specified by 
Dirichlet boundary conditions over some specified surface and there are many more points on such a surface 
than there are finite interpolation points of interest, it might seem that the result under consideration would 
be easy to show, but this is not the case and, moreover, there are various consequences of this result that 
are nontrivial (such as Theorem [5] in Appendix A and the fact that scalar point source basis functions are 
independent, which is otherwise hard to prove). 

Looking ahead, problems are encountered in the final stages of the proof when formal results dealing with 
outward continuation are required. Here outward continuation inside the interior of a sphere is the analog of 
downward continuation for the exterior of a sphere. Outward (or downward) continuation is well-known to 
be an ill-conditioned problem, at best. Consider, for example, the problem of outward continuation of / to 
r — 1, given that / is specified on some sphere of radius b, where f{b, 9, (p) — Ua{b, 9, (p) ioi Q < b < a < 1 
and where Ua{X) is specified in terms of a harmonic series for the same function introduced in the discussion 
after ([50|l . Although this leads to a solution that is formally harmonic since it is a linear combination of 
harmonic functions, one would not expect the partials of / for r > a to be bounded. The boundedness 
conditions on the class of admissible functions, 7i[o], rules out such examples. 

Starting with some / € 7i[o] specified over an appropriate interior boundary surface, if a formal solution 
to the outward continuation problem can be shown to exist, then provided that no solutions outside of 7i[o] 
are allowed to slip in there is a solution to the outward continuation problem, which is to say that the 
reconstructed harmonic extension of / is bounded and has bounded partials. Thus issues concerning ill- 
posedness versus ill-conditioning have implications that must be dealt with since an ill-conditioned solution 
processes may produce large but formally bounded results, while an ill-posed one, by definition, produces 
uncontrolled results. In this appendix an idealized theoretical perspective is taken so that it is considered to 
be sufficient to prove that something is possible, regardless of the associated implementation consequences 
(i.e., only the consequences of ill-posedness are relevant). (In the Section [5] actual implementation aspects 
are dealt with so the consequences of ill-conditioning must be, and thus are, considered.) The dividing line 
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here between ill-posedness and ill-conditioning is made here according to whether a constructive algorithm 
does or does not stay within the allowed family of admissible functions. 

In a general functional analysis setting, the problem of staying within some well behaved family of 
functions is well-known. The most common conceptual example encountered occurs when a sequence of 
smooth functions is considered that has as a step function as a limit, which is clearly outside the family 
of smooth functions itself. Here matters are a little more complicated since, besides 7i[o] , there are various 
relevant families of functions to consider (such as the family of harmonic functions, the family of smooth 
functions or the family of functions with bounded norm). 

It will turn out that, in the end, it is very hard to set up an independent constructive interpolation scheme 
that stays within 7i[o] when outward continuation is entertained. In spite of any difhculties encountered, 
observe that the stated theorems and interpolation results are obviously true since DIDACKS theory, in 
conjunction with the existing proof of point mass independence in the literature, immediately supplies a 
constructive harmonic interpolation algorithm that guarantees that end results are always in 7i[o]; however, 
the associated proofs arrived at by this route would not then be independent, which is the goal here. From 
the conclusions at the end of this appendix, in terms of characterizing outward (or downward continuation) 
as either an ill-posed or ill-conditioned phenomena, it would seem that it falls somewhere in between these 
two categories, being worse than what might normally be called ill-conditioned, but in some sense not being 
fully ill-posed either since a formal unique harmonic solution generally exists, even if it has unbounded 
partials. (It is this distinction that underlies the added boundedness requirements stated elsewhere here and 
in this sequence of articles.) 

To fix our goals, first consider an actual statement of the theorem of interest and for which an independent 
proof is sought: 

Theorem 8. 7i[o] uniformly pointwise independent: An f G 7i[o] always exists that assumes arbritrary 
(bounded) values at Nk < oo specified points, {Pk}k=i G ^[o],' i-S-, f{Pk) — Ck fork — I to Nk with arbritrary 
Ck 's. 

Observe that this theorem is about functions that are harmonic over iljo] — f2(o) U 9fi[o] and that a similar 
result obviously holds for functions that are harmonic over just fi(o)- Also, given any function that is 
harmonic over 51[o] it can be extended to a harmonic function over a slightly larger domain, so the issue of 
whether it is stated to include the boundary or not is not, in itself, especially germane (this is one of several 
reasons why it is assumed that admissible functions for the II norm are harmonic over i7[o] ^-i^d not simply 

^(0))- 

Obviously, Theorem [8] strongly suggests a proof by construction, but as just noted this does no mean 
that consideration must be given to practical numerical implementation issues. In particular no weight is 
given to, or needs to be given to, the consequences of ill-conditioning. 

Before proceeding, some additional mathematical apparatus and associated terminology is required. First, 
recall the usual concept of a star-shaped region, il, with respect to some specified interior point Xq: n is 
said to be a star-shaped region with respect to Xq if every straight-line segment connecting Xq and another 
arbritrary point of is constrained to lie wholly inside f2 (i.e., such line segments contain no points that are 
not in Cl). Here the primary interest is not in star-shaped regions themselves, but in the boundary surface of 
such regions; accordingly the term closed star surface will be used to denote the boundary of a star-shaped 
region. Moreover, since not only are the boundaries of star-shaped regions of interest needed, but also the 
hull-like counterparts of finite sets of points are needed, the required generalization will be labeled a star- 
surfaced set or a star-surface. Thus with this understanding, a general definition can be given that includes 
not only closed surfaces, but also finite sets of points. The criteria that a set of points E C K" (which can 
be either a finite or infinite) is star-surfaced with respect to some given point Xq is thus defined as follows: 
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Definition B.l The criteria for a region to be star-surfaced with respect to some given point Xq is that for 
all P € S, the line segment joining Xq and P must contain no other point of S. If S is a finite set or is only 
part of a closed surface it will be called an open star surface; alternatively, if S is a closed surface it will be 
called a closed star surface. 

Next, observe that this concept can be used to prove a theorem that will come in useful later (which is stated 
and proved for R" for n > 1, rather than R'^): 

Theorem 9. Let S C R", for n > 2 be a finite set of distinct points, then S is star-surfaced with respect to 
some point to or other, where to is constrained to be inside some given specified neighborhood o/R". 

Proof. Here it is obviously assumed that neither the neighborhood nor any of the points of E is at infinity. 
Let Nk be the number of points of S, then consider the set of extended straight lines L going through each 
pair of distinct points tk and t^' of E, which are indexed by the subscripts k and k': {Lkk'}- Let d{s, kk') 
be the (minimum) distance between a point s and Lk^t for all points of Lfe^' G {Lkk'}- Then consider some 
given fixed line L„ that is not in {L^k'} and let u £ L^, so that the function 

c?7r(u) := W d{u, kk') 

k,k'yik 

is zero when the point u and some points tk and tk' arc all common points of some line or other. (Here, 
while it is assumed that L„ ^ {Lkk'}, may be parallel to some element of {Lkk'} and in this case 
dT^{u) = const. 7^ may occur.) Since there are at most Nk{Nk — l)/2 zeros of d-,r{u) corresponding to the 
possible Nk(Nk — l)/2 intersections of the Lkk' 's with L„ there are at most a like number of coUinear points; 
hence, for any specified neighborhood it is obviously possible to find a point u d L^ where dj^iu) ^ and 
when to is taken to be such a point, S must be star-surfaced with respect to to. □ 

Given this result, the general strategy will be to first of all show that a closed star surface with respect to 
Xo can always be constructed from the given Nk points, where Xo can be taken to be a point to guaranteed 
to exist by Theorem [S] (Notice that it is not necessary here that Xo be in the convex hull of the original 
specified points.) Second, it will be shown that given such a surface, an interpolating function can be found 
that matches the prescribed point values on this surface; hence, this function also serves as an interpolating 
function for the given finite number of specified points. Third, it will be shown that this surface interpolating 
function can serve to specify a unique harmonic function over the interior of this region from the uniqueness 
of Dirichlet boundary conditions. Leaving, fourth and finally, the issue of performing an associated harmonic 
extension from this subrcgion to the whole of il[o]. Both of the first two steps here rely on the following 
lemma 

Lemma 1. Consider a set of distinct points {Uk}k=i that lie on a unit R^ sphere that has a center at Xq. If 
a set of bounded real (function) values o.fe specified that correspond to the locations {Uk}^Jli, then 

a smooth interpolating function F{X) (where \X\ = 1) can always be generated that matches the prescribed 
values: F{Uk) — Fk- 

Proof. The proof will be by construction. Let (r^, Ok, 4>k) = (1, 6k, (j^k) be the polar coordinates of Uk. Let 
the minimmn angular separation of the C/fc's be /3min- Also let the specified values of F be Fk := F{9k, <j>k)- 
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Further let Fq := Fk' be the minimum value of Fk, which is associated with the point Uk'- This value will be 
the default value of F{X) (i.e., the interpolating values at all points will take on this value unless otherwise 
specified). Notice that all that is really required is that around each interpolating point a smooth extrusion, 
which can be visualized as a local "function surface" deformation [where a function surface here is in the sense 
of ordinary calculus where z := Z{x, y) often occurs], be defined and added to Fq so that Fk will be matched. 
Next consider how this will be done for each point Uk in turn. To avoid confusion let k" be the index for the 
specific point under current consideration that is temporarily being singled out. First, temporarily erect a 
Cartesian coordinate system with a z-axis through the point Uk" . Reexpress the other 9^ and (jj^ 's in the new 
corresponding polar coordinates. Finally, for 9 < /3inin/2 let F{9k", <pk") = + {Fk ~ -Fo)(cos a{9)) + l)/2, 
where a{9) = tt{9/ f3ryiin), and F is otherwise as previously specified. After considering each point transcribe 
the results back to the original coordinates and repeat the procedure until it has been applied to all the 
points reindexed to k" . Although he interpolating function that results is clearly a smooth (not only C^, 
but even C°°) interpolating function, it is clearly not generally a practical one. □ 

The transition function here {i.e., H,^{X) = [cos(7rA") + l]/2, for X e M for 1 > X > 0} is known as a 
fifty-percent cosine taper. 

Notice that this lemma (and the following lemmas) can immediately be generalized to handle complex 
or vector valued functions by simply applying the lemma to each required component in turn. This lemma 
can also be used to show that a smooth interpolating surface always exists: 

Lemma 2. Given a finite set of distinct bounded points {Pk}k=i *^ r smooth interpolating surface exists 
that goes through these points; moreover, this interpolating surface is a star-shaped surface with respect to 
some point Xq so that the interior of this interpolating surface is a bounded star-shaped region. Moreover, 
Xq can be constrained to lie in any specified neighborhood. 

Proof. From Theorem|9]the set of points {Pk}k=i is star-surfaced with respect to some R'^ point Xq, which, in 
turn, is inside some specified neighborhood. Place the origin at Xq, and label all the coordinates accordingly 
so that, for example, 9k and cjjk will be defined as usual in this new coordinate system. These angles can 
then be used to define the unit sphere required by Lemma[TJ Uk = Pk/Pk, where Pk := \Pk\, as before. Next 
let Fk = Pk in Lemma[l] which guarantees an interpolating function F{9, 0) exists. Finally set r ~ F{9, 0), 
which yields a smooth interpolating function with all of the desired properties. □ 

Lemma [1] can be used again to show the following useful intermediate theorem: 

Theorem 10. Given a finite set of distinct bounded points {Pk}k=i *^ '^'^'^ associated (function) 

values at these points {Fk}^^^, a smooth interpolating surface S exists and a smooth interpolating function 
F{X) also exists such that Pk G S and F{P^,) = Fk for all k; moreover, F{X) is smooth for all X ^ Yi and 
S is star-shaped surface with respect to some point Xq, that lies inside a given neighborhood, which can be 
taken to be one that also contains the given coordinate origin. 
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Proof. Lemma [5] immediately yields an interpolating surface E with all of the desired properties, but it 
remains to be shown that an appropriate F(X) exists. First reposition the origin over Xq, then consider a 
function G = G{9, (p) defined on the miit sphere such that G{9k, (t>k) — Fk- Lemma [T] guarantees that such a 
function exists. Finally, since S is a star-shaped surface F{r, 9, (f) — G{9, 4>) is a well defined function that 
has all of the desired properties. This function can then be transcribed back to the original coordinates. □ 

Taking stock of where we stand. Theorem [TUl shows immediately that given a set of points {Pk}^^i in 
ri(o) ^ smooth interpolating surface E and function F{X) exist for X e E. From the uniqueness of Dirichlet 
boundary conditions, it follows that a harmonic function f{X) exists that is defined throughout the whole 
interior of E and which takes on the values f{X) = F{X) when X e E. What remains to be examined is 
the possibility that the specified boundary values on E also can be harmonically extended outward to cover 
the region between E and dQ[o]- If this can be done, then it will show that a smooth closed interpolating 
surface E exists with smooth specified boundary conditions F{X) and that it can be uniquely extended to 
yield a function that is harmonic throughout all of r2[o]. Towards that end the following lemma is useful: 

Lemma 3. For f G 7i[o], */ /(^) = throughout some sphere of diameter less than one that is centered 
over the origin o/il^o]? then f{X) = throughout Q[q]. 

Proof. Following the notation of Jackson [11], let Yim{9, (f>) denote the (complex) spherical harmonics spec- 
ified by 

where Pl"{x) are the associated Legendre functions. For complex constants aim the Dirichlet boundary 
value problem for ?i[Q] in fijo] has the solution 
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and so all / € 7i[o] also have this form. (Although the aim and Yim are complex valued they can obviously 
be used to represent the real boundary value problem here also.) 

Let a denote the radius of the sphere that is centered over the origin, then evaluating (jB-2[) at r = a and 
using the ortho- normality conditions for the Yi ,„ , it follows immediately that ai m = for all allowed values 
of I and m. Hence from the form (jB-21) itself / = throughout f2[o]- □ 

This result can easily generalized: 

Lemma 4. For f G 7i[o], if f{^) = throughout some sphere of finite diameter enclosed wholly inside f2[o]; 
then f{X) = throughout all o/iljo]. 

Proof. Again let a be denote the radius of the sphere in question. A sphere of radius a/2 can be situated 
wholly inside this given sphere as close to the origin as possible. Consider a spherical harmonic expansion 
like (jB-21) . except with spherical coordinates that have an origin at the center of this new sphere of radius 
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a/2. Evaluate this new expansion in exactly the same way as in Lemma [3] This harmonic expansion can 
then be used to show that / = thorough a sphere of radius (3/2)a centered over the sphere of radius a/2. 
Next again place a sphere of radius a/2 wholly inside this new sphere as close to the origin as possible and 
repeat the same argument. After a finite number of iterations, clearly the origin will eventually be enclosed 
by this sequence of new spheres. As soon as this happens the required result immediately follows from 
Lemma [S] □ 

Lemma |4] can be used to show that a harmonic function specified over an interpolating surface that is 
inside 51(0) has a unique harmonic extension throughout the whole of ri[o]- First Lemma |4] implies: 

Lemma 5. For f e 7i[o], if f{X) is specified over some sphere of finite diameter enclosed wholly inside 
f2(o), then it has a unique harmonic extension throughout the whole o/fi^o]- 

Proof. Suppose to the contrary that this is not true and let S denote the sphere in question, then / and 
g e H[o] exist such that g{P) = f{P) for all P e S, but g{X) ^ f{X) for some X £ fl[o]- Let x{X) = 
f{X) — g{X), then from the fact that x{P) = for all P e E and the uniqueness of Dirichlet boundary 
conditions Lemma [3] can be invoked to conclude that xi^) — throughout all of fi^g] and a proof by 
contradiction follows. □ 

Since the interpolation surfaces constructed earlier all enclose a sphere that, in turn, is inside ri(o), it 
might seem that the desired general result immediately follows, but it is at this stage that the proof falls 
apart because of a hidden assumption that was hinted at in the beginning of this appendix. Before pointing 
it out and drawing the appropriate conclusions, it is useful to summarize the results so far. 



Summary and Conclusions 



In summary, for a specified interpolation procedure, it has been shown that a set of interior points can be 
considered to be star-shaped and that an interpolating function to the given /(Pfe) values can be constructed, 
which is unique over this surface. The uniqueness of Dirichlet boundary conditions thus imply that the values 
of / are uniquely determined throughout the interior of this star-shaped region by this interpolation scheme. 
Clearly, the values at each of the points /(Pfe) can be assigned differently and independently from each other 
by the very nature of Dirichlet boundary value problem itself. Moreover, a sphere of finite radius can be 
wholly contained within this star-shaped region and / can be uniquely specified over it. In turn, by Lemma[5] 
this / can then be harmonically extended to the whole of ^[o]. Thus our initial goal — an independent proof 
of Theorem |5] — would seem to follow; however, there is a hidden assumption that has been made. Namely 
that a smooth function (i.e., C°°) specified over a smooth boundary specifies not only a harmonic function, 
but one that is bounded and has bounded partials (i.e., a harmonic extension in 7i[o]-) To see the distinction 
here consider let E be a sphere of radius 6=1/4 and let ai = 1.01 and a2 — .99 where ai and 02 are specified 
a parameter values for the function Ua discussed at the first of this appendix. Let Uai and denote these 
two functions with a = ai and a = a2 in UaiX), respectively. Then it is clear that Ua^iP) ~ Ua2{P) for 
all P S E and that both functions are C°° on S. Furthermore, the values on E have a unique harmonic 
extension to all of r2[o] (namely Uai and Ua2)] however, [7^2 G W[o]i but Uai </ W[o]- Finally, it seems to 
be hard to engineer a constructive interpolating function, with an associated interpolating surface, that 
is simultaneously guaranteed to match the specified values and be in 7i[o] — although this leads to some 
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interesting research problems in itself. [For example, what are the N pointwise independence properties of 
truncated spherical harmonic series (including geographic factors) and how can they be proved?] 

These general sort of examples have been discussed elsewhere within a downward continuation context 
(see, for example, [17 ) and they seem to indicate that some sort of direct regularization is called for when 
implementing harmonic extrapolation. In particular, although Theorem [T] guarantees that a weighted en- 
ergy minimization principle underlies DIDACKS interpolation, which in itself is a very useful criteria that 
guarantees an answer in 7i[o] , if data errors are present or if the source locations are inappropriate then the 
implied harmonic extrapolations may be, at best, useless. [Here even in the presence of exact data, at a 
minimum Nyquist frequency and other related digital signal processing issues (such as aliasing) probably 
should always be considered.] 

Nevertheless, several useful intermediate lemmas and theorems have been presented in this appendix. 
Furthermore, as previously pointed out the Kelvin transformation can be used to obtain corresponding 
results for the exterior of a unit sphere and many of the lemmas and theorems clearly apply to regions of 
a more general shape. Also observe that these interpolating results also apply to analogous situations in 
M^, where the lines-of-proof followed are even easier to follow. Many of the results here can also be readily 
generalized to different shaped regions in M" for n > 3; however, some technical issues must be dealt with 
to prove corresponding theorems in these general settings. 
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